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INTARES AND OUTELESS POR LOW- HEAD HY DROPOWER = 


power generation; Low head; Piers; Standardization; 


development is presented. In addition, the current flow passage design practices are 
sumnrarized. Several package units are available in the range from 100 hp to 6,700 hp © 
(75 kW to 5,000 kW). These predesigned units reduce equipment costs through 
standardized manfacturing techniques. Several types of package units are considered. _ 
« Manufacturer’s standard flow passage dimensions are given for bulb and rim-generator — 
_ turbines. These standard dimensions can be used for initial layout of a low-head 
hydroelectric development. Standard flow passage designs are feasible for similar sites; — 


site conditions vary. Simplified intake and draft tube may less costly 


' _ REFERENCE: Pugh, Clifford A. (Hydr. Engr., Dept. of the Interior, Water and Power © 
Resources Service, Engrg. and Research Center, Denver, Colo. » “Intakes and Outlets 
H for Low-Head Hydropower,” 
| HY9, Proc. Paper 16494, September, 1981, pp. 1029- 1045 


‘KEY WORDS: Algorithms; programs; Equations; 
Pipe flow; Pumps; Regulators; Reservoirs; —— distribution ' 
ABSTRACT: Unknown pipe be helpful in the design and operation 
| of water distribution networks. By using unknown pipe characteristics it becomes easier _ 
‘ to select appropriate pipes, to choose boosters, and to regulate the flow through pumps — 
using regulating valves. An algorithm is presented for automatically constructing and 
= equations describing the behavior of water distribution networks when some 


pipe characteristics are unknown. It is based on loop equations and is therefore 


of the Division, 107, No. HY9, Proc. 


16492 PARAMETRIC STUDY OF FLOOD WAVE PROPAGATION 

KEY WORDS: (waterways); Finite difference method; Floods; 

Flood waves; ‘Hydrographs; Open channel flow; Parametric hydrology; _ 


ABSTRACT: A parametric of the flood wave propagation. probles | is 
based on numerical solution of the nondimensionalized unsteady flow equations of — 
ig channels. The propagation of a sinusoidal flood wave in a prismatic channel : 


studied for uniform initial flow. The governing parameters (initial uniform flow Froude 
1 number, wave amplitude, wave duration, channel width parameter and side slope) are 
z varied over a wide range. In all, 49 cases are studied. Effects of these governing 
_ parameters on the subsidence of stage and discharge and the speed of the wave peak 
- g are described in detail. The relative wave amplitude is found to vary linearly with Fo 
ig the initial uniform flow froude number, for lower Fy values. Wave duration has a very 
aa sa effect on _ Subsidence a greater subsidence at lower wave duration 4 
REFERENCE: Sridharan, K. Prof., Civ. Engrg. Dept., Indian of Science, 
= Bangalore, India 560012), and Kumar, M. S. Mohan, “Parametric Study of Flood 
Wave Propagation,” Journal of the Hydraulics Division, HY9, 
Proc. Paper 16492, September, 1981, pp. 1061-1076 


one A state-of-the-art review on standardization in low-head hydroelectric 


however, structural and geological considerations may cause design variations when a 5 
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16506 DEMAND FORECASTING IN WATER SUPPLY NETWORKS 
KEY WORDS: Computerized scheduling; Computer programs; Consumer 
demand forecasting; Demand (economics); Filtration; Optimization; Rie 
Predictions; Time factors; Water demand; ‘Water Water Water supply 
forecasting 
ABSTRACT: The problem of water demands over a 24 hour 
period is formulated within the framework of an on-line computer based pump _ 
scheduling scheme for optimum operational control of a water supply network. Two — 
_ well known forecasting techniques are considered and compared, in terms of accuracy : 
and computational using data from a network in in United — 


REFERENCE: Perry, Patrick F. (Lect., Engrg. Dept., “Imperial Coll., 


SW7 2AZ, England), “Demand Forecasting in Water Supply Networks,” Journal of the 

Hydraulics Division, ASCE, Vol. 107, No. HY9, Proc. Paper 16506, September, 1981 


vs 
16505 SEDIMENT-CONTAMINANT TRANSPORT MODEL mes 
KEY wW ORDS: Contamination; Estuaries; Finite element method; — 
materials; James River; Models; Rivers; Sediment; Sediment 


_ ABSTRACT: The unsteady, two-dimensional model, FETRA, was developed to 
_ simulate both sediment and contaminant transport in rivers and esturies. The model 
_ consists of three submodels which, when used jointly accurately depict the interaction : 
and migration of sediment, dissolved contaminants, and particulate contaminants. 
_ FETRA solves the migration (transport, deposition and resuspension) of cohesive and 
_ Roncohesive sediments and particulate contaminants of three sediment-size ag 
FETRA was applied to the James River estuary to simulate the migration of river 
sediments and the pesticide, Kepone. Computed results produced by the model 1 = 
_ very similar to field-measured data, meaning that FETRA is reasonably accurate. ay 
t REFERENCE: Onishi, Yasuo (Staff Engr., Pacific Northwest Lab., Battelle ls 
Institute for the Dept. of Energy, P.O. Box 999, Richland, Wash.), . : 
Contaminant Transport Model,” Journal of the Hydraulics Division, ASCE, Vol. 107, 
No. HY9, Proc. Paper 16505, September, 1981 1981, 
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NTAK S AND OUTLETS FOR Low-Heap 


By Clifford A. Pugh," A. M. ASCE 
as or both, resut ta hese ope: Nenedity, “The 
rth package | the packages aca, 

_ The need for additional energy production in the United States ial socal resulted 
in increased interest in the development of remaining hydropower resources. 
Much of the undeveloped capacity is in the low-head range (less than about 
_ The main obstacle to development of low-head sites has been economics. 
Te. cost per | installed kw (kilowatt) is still higher than for fossil fuel plants 

in most cases. However, with the cost of fuels continuing to rise, the low-head © 
“hydropower alternative is becoming more favorable. 
_ This paper outlines present design practices regarding flow oeens tis 
to determine if standardization or design changes are possible to reduce construc- _ 
tion costs. In low-head plants, head losses could reduce the net effective head — 
significantly. Therefore, it would be desirable to streamline the flow passages 
as much as possible. At the same time, economics dictate that the flow passage 
should be simple and small. With these conflicting interests in mind, the present 
design methods were examined to determine possible design changes to simplify — 
present designs without introducing significant additional losses. 7a ia 
The design of penstock entrances for high-head dams is one example that 


— 


illustrates the value of taking a close look at. conventional design ‘practices. 2 


It was conservatively estimated that $13,000,000 was saved in construction costs 
_ on the penstock entrances for the Third Powerplant at Grand Coulee, by reducing 
the size of the bellmouth entrances (18,19). Penstock entrances have historically — 


been designed using the same criteria used in designing high-velocity conduits, — 
_ whereas the velocity | in the -penstocks is much lower. 


This state-of-the-art survey is based on information obtained in a literature 
_ search and on comments made by manufacturers and consultants in low- head 


“Presented at the August 8- ASCE Divisions Cae, 
_ San Francisco, Calif. 
"Hydr. Engr., Dept. of t the Bureau Reclamation, Engrg. and Research arch Center, 

_ Note.—Discussion open until February a 1982. To extend the closing date one month, 

a written request must be filed with the Manager of Technical and Professional Publications, 

_ ASCE. Manuscript was submitted for review for possible publication on October 17, 
1980. This paper is part of the Journal of the Hydraulics Division, Proceedings of the 
_ American Society of Civil Engineers, ©ASCE, Vol. 107, No. HY9, September, 1981. 
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_ In the head range between 10 ft and 148 ft (3 m and 45 m), turbines with 7 +4 
_ propeller-type runners are typically used (12). If the Tunners are adjustable, — 
the turbines are called ‘Kaplan turbines.”’ These turbines are usually om 

with a vertical shaft, a spiral case, and an elbow- -type draft tube (Fig. 1). 

ss A large percentage of present low-head turbines are this type. p at 
heads less than 66 ft (20 m), axial flow turbines have proved to be more econom- — 
‘ical. This is concerned d with the the low- head omg, than 66 


a ft (20 m). 


» 


FIG. 1, Kaplan Turbine 


Pooper == 


Turbines in which the water is conducted to the distributor coaxially with ; 
_ the shaft are called ‘‘axial flow turbines’’ or sometimes ‘“‘tubular turbines.’ 
7 To avoid confusion, the term “‘axial flow y turbine’’ will be used in this paper. 
Axial flow | turbines also use propeller- type runners. In s some | cases, adjustable 
(Kaplan) runners are used. The three primary types of axial flow units are 
_ Shown in Fig. 2. 


| 
err 
| 


= few package-type watied: and generator u units are available in the low- head 
range. These predesigned units reduce equipment costs by eliminating the need 
for site- specific engineering and by using standardized manufacturing techniques. , 
However, the standard flow passage shapes | result i in a loss in turbine efficiency. j 
Butterfly valves (often used in the intake of package units) increase losses 
_ and cause an uneven flow distribution at the runner resulting in reduced operating — 
efficiency. Standard-draft tube shapes also lead to a loss in efficiency due to 
' 7 sharp corners. Package units have either fixed runner blades or fixed wicket 
gate positions or both, resulting in less operational flexibility. The economics Ca 
4 of each potential site should be evaluated to determine if the savings obtained z 


- with package units outweigh the losses. The following package units are pt 


¥ > Chalmers Tube Turbine ‘Units. - Ten standardized packaged designs ¢ are 


_ 49 ft (1S m). Flow is controlled with a butterfly valve in the intake. The es 
gates are fixed and the runners are either fixed or adjustable. ais hl 
Karlstads Mekaniska Werkstad (KMW) Miniturbines.—These units are 
for horizontal or vertical installation for flows from 35 /s—530 ft’ /s (1 m*/s-15 
= m /s) and heads from 13 ft-82 ft (4 m-25 m). Outputs range from 134 hp-2, 413 
h hp (100 kW-1,800 kW). The turbines have fixed guide vanes and fixed runner 
blades, although the position of the latter can be changed when the turbine _ 
_ is stationary. Flow is controlled with a butterfly valve. Fig. a aa a standard 


- Fuji Package- type Bulb Turbine and Generator.—Available in 19 models 
covering a range of net heads from 16.4 ft-51 ft (5 m-18 m) and outputs from a 
- 402 hp-5,362 hp (300 kW-4,000 kW). The runner is of the fixed blade type Jf 
and the flow is controlled with movable wicket gates. A sectional drawing of 
the standard bulb te turbine and the generator is shown on Fig. 4. 
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wae. 
Neyrpic Standardized Right Angle Drive Unit. —Units vary in output from 
“100 kW-1,500 kW and cover a head range of 6 ft-60 ft (2 m- -18 m) with discharges _ pee 
of 35 ft? /s-750 ft*/s (1 m°*/s-23 m’/s). Flow is controlled with a downstream 
control gate at the end of the draft tube. Seven standard units are available 
with runner diameters from 18 to 71 inches (450 mm- 1,80 800 mm). m2 8 shows 


a cross section through a right angle d drive tt turbine. 
‘cross flow”” * turbine assemblies are available in a range of units 
_ which have the ability to ‘cope with variations in head and discharge on small © 


DISCHARGE RING 
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3 BULB FOR TURBINE 
WICKET GATE 
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dams. The cross flow turbine is a radial impulse-t -type ——, The water is 
forced through a guide vane system and the blades of the cylindrical runners. a 
The water then passes through the runners again to the outlet. Fig. 6 is a 
cross section through a cross flow turbine . The adjustable guide vane allows. 
power generation from 16%-100% of design flow. Although the efficiency of — 7 
this unit is not as high as for axial flow turbines, it may have widespread — 
- application to small sites where limited storage is available and flow and head — 


| 
| 
| 
; 


7 A few manufacturers provide drawings with standard flow passage dimensions. 
- These dimensions were determined as a result of model tests and their operational : 


experience | over a period of years. _ Fig. 7 shows a - cross 


section and plan of 


ADMISSION 


bulb turbine according to Escher The are given in terms 


of the runner diameter D-1. Fig. 8 shows Escher Wyss’ standard dimensions re 
for the ‘‘Straflo’’ rim-generator unit. Fig. $ 9 shows Fuji’s s standard flow passage — 


= 7 By) res 
= ~ VERTICAL ADMISSION 


At this time there are no universally accepted flow passage shapes. Studies . 
on optimization of flow passage shapes done by manufacturers are considered 
poopelotary information. Therefore, the effect of intake and 
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1G. 8. —Standard Rim Flow Wyss) 
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i“ FIG. 9. —Stenderd Bulb-Turbine Flow Passage Dimensions (Fuji Electric Co. 
widely among manufacturers. The standard flow passage designs according to 
Ss and Escher Wyss for bulb turbines illustrate the large variation in design _ 


the civil works. 


— 
| 
flow passage shape to hydraulic efficiency, a study could be done to evaluate - 
: ess of simplifying the intake shape to reduce the cost of! 


INTAKES AND OUTLETS 
Stanoano DesiGns ‘FOR Line OF 
9 _ An alternative to standardization of flow passages in general would be to. 
develop a ‘standardized design for a series of powerplants to be installed at . 
similar locations. A feasibility study was performed by Motor-Columbus Consult- 
a ing Engineers to investigate harnessing of the Rhine River upstream of Lake 
_ Constance. The study investigated using 46 standard bulb units in 16 stations. 
Fig. 10 is a general layout developed for the project. This concept may be 
ee to development of low-head power on American rivers with several 


_ possible low-head sites, such as the Ohio or Mississippi Rivers where numerous 
navigation dams without power ——— currently exist. A standard design 


Stations in Sw WEDEN els 
_ The Swedish Power Association has initiated a program to overhaul and replace 
= discontinued small power stations 134 hp-2,011 hp (100 kW-1,500 kW) in Sweden > 


aw 


FIG. 10. —Powerplants on Rhine River from Domat (EMS) to Flasch (Motor Columbus 
The program includes hundreds of stations having a total installed capacity 
= about 3. 15 hp > x 10° hp (235, MW). few simplified automated axial flow 


; the 1975-77 “period to study operation 1 and economics s of the units selected. iarins 
In most of the small installations in Sweden it is possible to run the power 
sation intermittently. This makes it possible to design the turbine with fixed — 
_ blades and fixed guide vanes and omit the regulator. The turbines are 
mun as the reservoir is drained and shut off as the reservoir fills up again. 

oy If the installations cannot be run intermittently, they can be run with regulation 
for constant water level upstream using units with movable runner blades. Turbine 


The pilot sta stations indicated that i it was difficult to the costs down 


| 
: _ has made economic aid available (up to 35% of the cost) for the least profitable ; 
7 projects. Economic aid was considered to be justified since hydropower reduces 


“dependence on on oil and is acceptable 
form of energy that would otherwise be wasted. __ 


Most designers and consultants contacted feel that standardization of the 
inlet and outlet sections for low-head hydropower installations will not always — 
be possible. The shape of the water passage is influenced by the structural 
= system requirements and the geology of the site. Even units designed 

by the same manufacturer vary from project to project. 


7 - This section summarizes material in technical literature, and comments made 
by manufacturers and designers concerning flow passage design. Many points 
Bes in the references are described; however, the Teferences should be 
referred to if more detail is required. 
Forebay ._—The design of the entrance channel in the forebay should provide 
a uniform flow distribution to avoid the tendency for vortices and to minimize 
~ trashrack losses. If the geometry of the approach channel is in question, a 
model study can be done to assure uniform approach flow conditions. Flow 
velocities and surges are of particular concern on navigable rivers. -Run- of-river 
_low- head dams are typically designed with the powerhouse on one side of the _ 
_ river and an overflow weir-type dam on the other side. In some cases an entrance 
_ channel to the powerhouse is excavated in the bank. This may cause problems 
ol of swirling flows in the approach channels and cross flows in the exit channel. | 
_ Intake—Vortex Formation.—Intakes should be designed with a sufficient 
; submergence to avoid vortices. _ Air- -entraining vortices may decrease turbine 
> efficiency and pull floating debris into the turbine. Gordon (7) describes the 
development of design criteria to avoid vortices at low-head intakes based on > 


a study of 29 existing hydroelectric intakes, 


‘The intakes studied have the general configuration shown in Fig. 11. 


factors which : appear to affect the formation of a vortex. are: (1) Geometry — — 
of the approach flow; (2) velocity at the intake; (3) the size of the intake, a 


and (4) submergence. This article concentrates on the effects of velocity, intake | 


size, and submergence on vortex formation. Gordon developed the following © 
The ‘effects of is approach flow were not evaluated in this study. Ho 
; Montreal Engineering Company of Canada uses a coefficient C of 0.3 for intakes _ 
with symmetrical flow and 0.4 for lateral approach flow (7) to determine minimum Fad 
e allowable submergence. These coefficients correspond to the lower and upper 
limits, respectively, of the shaded area on Fig. 11. 
If the available submergence is not adequate, a model study can can be done 
‘to study vortex formation. However, accurate prediction of vortex formation 
using a scale model is difficult. The effect of surface tension and air entrainment P 
are important and are difficult to study in a Froude scale model. Durgin (6) 
es a general n method for i aaa potential scale effects of free surface 
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vortices, so that a projection to prototype operating conditions can be made. only 
_ Zeigler (27) investigated the use of rafts placed under the water surface and : 


other devices to prevent formation of air-entraining vortices at Grand Coulee 


Third Powerplant. The hydraulic mc model study included vortex tests on the effect 
of trashrack structures, upstream channel geometry, intake modifications, 4 


_ Bisaz (3) reports on the use of the flow ‘“‘injector shaft.’’ This device has 

the same effect as a raft. Because of an increase of surface current, vortices 

: are suppressed. The injector was developed for bulb turbines on run-of-river 
plants to increase flow velocity on the surface to promote movement of floating _ 
_ debris toward the trashrack on the front of the powerhouse. The configuration — 


of a typical flow injector is shown on Fig. 10. ides T 


¢g Pier Design.—In typical run-of-river plants a flow separating pier divides the 
powerhouse and spillway. Piers also divide the intakes where there are multiple 
units. The size and shape of the x pier should be carefully designed to avoid — 


= areas or flow disturbances which may cause excessive head losses 


é CONF AN INTAKE 
with vortex problems 
intakes with no vortices 
Recommended min submerge 


RECOMMENDED MINIMUM SUBMERGENCE 


Corps of Engineers (22) contains design information useful in designing piers © 


4 Rouve (20) outlines German practice in designing the size and shape of flow > 
separating piers” between a powerhouse and weir. The pier design procedure 
takes into account the location of the powerplant with respect to bends in o 

= the river, alinement with the weir, and the entrance channel to the powerplant — 
into the riverbank. Fig. shows the general shape of the flow 


the following | formula defines the width of the (B): 


Ne? 
The coefficient C varies between 0.7 and 1.4 depending on the mentioned 


| 
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; omen Ref. 20 should be referred to for det details of this design method. 
The HDC sheets 111-5, 111-6, and 122-2 (22) describe pier effects for gated im 
“overflow spillways. Abutment effects are described in sheet 111-3/1 and 111-3/2. _ 
‘ When spillways are operated with adjacent bays closed, the piers adjacent to i 
_ the closed bays produce abutment-type effects. Although these pier design criteria — 
are for ‘Spillways and not for submerged intakes — do show the relationship | 


13 ‘Shows five pier 


low heads. Types 2, 3, and 3A were recommended for general use with high 
7 heads. Type 3A had the most desirable flow contraction characteristics.  =— 
g Trashrack Design. —Trashracks are typically used to Prevent large debris from 


rectangular bar trashracks. The article. evaluates previous ; work on trashrack 
and baffle losses and introduces additional data to supplement design information. , _ 

Equations are developed to calculate the head loss as a function of the approach 

and clear between the The head loss through 


4 


= 
-12—Flow Separating Pier betwee FIG. 13.—Pier Nose Shapes (Type 3A 


> 


rd Dimensions in | Parentheses) 


the trashrack was found to be a minimum for bars with depth over thickness - . 
“@ ratios s of about 3. 0 for any solidarity (S) between 0. 25 and 0. 33, —— 


possibility of flow induced trashrack vibrations should also be 
in the design. Trashracks spanning large intakes or subjected to high velocities — 
are more likely to experience vibration problems. Vigander (24) describes — 
_trashrack vibration studies on the Raccoon Mountain Pumped-Storage Project. 
ott was concluded from the tests that damping by either rubber pads or hydraulic 
baffle plates would effectively reduce flow induced vibrations on the Raccoon 
ae Shape.—Intake shapes vary with the type of | turbine 1 used. Intakes 
are either circular or have a rectangular to circular transition. The face of the 
_ intake is either vertical or inclined. The hydraulic design of the intake is based 
_ largely on engineering judgment; the classical bellmouth criteria are used as 
@ Starting point. The entrance end of the bellmouth curve is then truncated — 7 
2 to obtain a short intake profile which achieves required contraction in a bongeh 7 


nose shapes. Pier types 1 and 4 are the least  &§ 
__of negative nressures Owever nevative pressures may not he a concern for 
| 
ay 
$i 
FIG 
| 


shape on losses is not Studies by Peterka (18) 
* Rhone (19) indicate, however, that the size of the intake could one -vdl 

te reduced without introducing excessive additional hydraulic losses, since flow 


separation and resulting energy losses are less of a problem at low velociion. 
_ The ideal bellmouth entrance design was developed for high-velocity conduits. j 


Since the velocity for low-head turbines is much lower, it is possible that savings 
could be achieved in trashracks, bulkheads, entrance gates, and other areas” 


itt might also be possible t to > simplify the framework by using { ~~ ‘Surfaces 


changes require further ‘study to determine. if reducing entrance sizes 
intake shapes are worth pursuing in low-head powerplants. 


KMAKNAN TT 
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—g Fish Passage Allowances.— ntake \ velocities ‘may be limit limited if if fish ilies 

d facilities or screening devices are included in the design. This limitation may 
control the intake design. Optimum water velocities will depend on the type 
a fish and on the type of facilities used. Bell (2) gives fishway structure ire design 
criteria and other factors related to design of fish passage f facilities. — ‘dhe. Ca 

w Downstream passage of fish through low-head axial flow turbines may not 
be a serious problem since the openings between the wicket gates and the 
runners are generally large, the velocities are low, and the flow passages are > 
‘Straight. Bell (2) states that “‘turbines of modern design generally have e fish 
4 Runner Diameter and Turbine Setting.—The turbine setting is determined by 
th 


e allowable specific speed which i is susnally controlled by the cavitation potential. . 


‘ied 
q 


an economic it is desirable t to select the speed 

possible since this results in a smaller, lighter generator. The cavitation potential | 
of an axial flow turbine is generally lower than a Kaplan unit. This allows 
~ axial turbine to have a higher setting or a higher speed (12). Other factors — 
such as smooth operating range and submergence also enter into the design. | 

Heinemann (9) used data from 130-power stations to develop a design aid io 

e determine runner diameter and the powerhouse size. The upper chart in Fig. p 

: 14 makes it possible to determine the capacity in kilowatts. In the lower diagram — 

i the runner diameter can be determined. The head, H, is marked on the lateral 
_ scales and a line is drawn connecting these points. The runner diameter can 
then be determined for a certain discharge (Q) or | speed (n). The range of | 

_ applicability is bounded by a region of uneconomic operation and the region 

of cavitation. The assumptions were — gi comparing existing al 
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_ Sutherland (21) also gives charts for finding the approximate runner diameter - 
and rules of thumb enabling the approximate dimensions for an entire unit 
7 to be found. An empirical relation: between diameter (D) and the ratio of f power 


to head lo head (in h horsepower / head) i iets 
tf. metric units. For power in terms of kilowatts of turbine output the formula 
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Dadu (4) presents a developing 2 a for determining the 
Setting of a hydraulic turbine relative to the tailwater level. The turbine 
manufacturer must supply data concerning cavitation performance to develop ; 
_ the nomograph. This method is useful only for the installation under consideration. | 


= Draft Tube.—The hydraulic design of the draft tube is customarily ro , 
7 
> 


by the turbine manufacturer, since it is considered an integral part of the turbine | 
in determining the turbine performance. However, there is potential for decreasing 
the length of draft tubes, thus decreasing the cost of the civil works. According _ 
- to Kline (11), the most efficient conical draft tube would be about 10-runner 
_ diameters long with a total included angle divergence of about 7°. Economics 
_ _ usually limits the length of the draft tube to 4.5-runner diameters—S-runner 
a diameters with the total angle of diversion of about 13°-15°. A study done 
on the feasibility of low-head hydroelectric generation in 19 1969 by Mercer (13) 

= concluded that, ‘‘The section of the civil works structure where improvement _ 
would be most significant is the draft tube.’” The draft tube accounts for about 

_ 30% of the cost of the civil works in a low- head structure. There are several — 
possible methods presented by Mercer for preventing boundary 


thus shortening the draft tube (Fig. 15), 


ad 


adi 


4-ExIT VELOCITY ROTATION MOSSES. 


16.— —Model Bulb Turbine Balances 


‘The study concluded that a need existed for research in the design of draft * 


~ tubes to apply modern principles of boundary layer control to reduce the length — 
g draft tube diffusers. Research in this area has not been done on hydraulic — 
turbines. The concept of boundary layer control has been bypassed in favor 
: of simple conical diffuser designs because of the potential risks of complications —— 
a involved with a new concept. Shortening the draft tube in this manner would 
probably require a successful demonstration project before it would be accepted. | 
(15) reports on research to shorten diffusers for wind turbines using» 
boundary layer control. This research indicates that there is a promise of reducing © 
the length of diffusers by an order of magnitude while maintaining efficiency. vi 
jee (25) deals with the question of whether addition of swirl to fluid, 4 
fa flowing axially, in a diffuser would improve its performance. Model 
studies indicated that the coefficient of performance of conical diffusers is 
improved by an addition of a swirl velocity component corresponding to a forced — 
vortex. This paper concluded that the optimum swirl angle is about equal to 
the diffuser’s total conical angle. Improvements are not likely to extend to 
diffusers with total conical angles greater thar 30°20 
Varlamov ‘feports on a study that investigated interaction 


| 
| 
| 

4 | 
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2 turbine performance. Results of model tests done by manufacturers are not 


as flow moved ‘through — various s elements 0 of met flow f passage ina bulb 
turbine. Flow parameters were measured in various turbine operating modes, 
characterized by blade angles, 6 = 0°, 5°, and 20° and rotating speeds N = 
110 ) r/min, 140 ) r/min, and ‘180 r/min \ with various $ settings of the wicket t gates. 


i for two turbine operating modes. It should be noted that the selection of design 
data should be based on minimum total losses in the turbine as opposed to 
minimum losses in each section of the waterway. Fig. 16 shows the dimensions | 


_ are presented ina chart for several operating modes. Fig. 16 gives energy ce 


Alestig (1), de Siervo (5), (10), Mikhailov (14), and Worster (26) provide 
additional information on draft tube design. Gubin (8) describes developments 
in design and construction of draft tubes in the Soviet Union and western countries. 
Several consultants indicate that the section of the draft tube immediately 

_ downstream from the runner is the most important. The design of this section y 

_ is critical to the efficiency of the draft tube. Therefore, draft tube gates, and 
other obstructions should be at least 2-runner diameters downstream from the =. 


‘from the oun tube at low ealtiousien 3 may cause sneges and turbulence and the 
er level. These 
‘tailwater. fluctuations can c cause | power swings in the turbine The 
channel and the tailrace should be designed to minimize these problems. However, | 
tailrace design guidelines are sketchy, consultants recommend that a steep rise 


in the channel bottom be avoided, but the allowable rate of rise is not defined. a 


‘Complete of flow passages for low-head 
ments is not feasible. Structural and geological considerations may cause design 
_ variation from site to site . Allowances for fish passage may control the design 
in some cases. Package units are presently available for small plants. These : 
: units have standard flow passages included. The low costs resulting from — 
standardization are achieved at the expense of some loss in efficiency. However, 
the cost savings may be significant especially if several installations are considered 


Standardized designs are for similar site characteristics. Where geologic 


fora a series s of installations, , thus reducing the « costs s associated with | site- -specific - | 


_ The hydraulic design of the draft tube is customarily determined by the 
‘manufacturer since it is considered an integral pait of the turbine in determining — 


- publicly available. However, modern techniques of boundary layer control may = | 
_ be useful in shortening draft tubes, thus considerably reducing construction — 


= 
— 

= @€=€=—S- The outlet from the draft tube is generally designed to have a ~ Fy 

: = 

&g 
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The ‘possibility of reducing draft tube lengths through boundary layer 


has not been investigated in hydraulic turbines. 
Intake shape is typically determined using design criteria for high welachi 

‘conduit entrances. It is possible that intake shapes may be reduced in size 

or simplified without a significant loss in efficiency. Velocities are generally 


low in turbine intakes and the resulting | energy losses are also low. Simple 
A flow passage — can be used in these low velocity areas without increasing 
5 The work upon which this paper was based was sponsored by the Department 
of Energy under agreement EG-77-A-36-1024 with the Water and Power Resources" 
Service. The following companies provided some of the information used: Harza e 
Engineering Company, International Engineering Company, Allis-Chalmers, 
Motor Columbus Consulting Engineers, Inc. . F. W.E. Stapenhorst, Inc., Fuji 7 
Electric Company, Ltd., Escher Wyss, Public Utility District No. | of 4 
‘Chelan County, Washington, Karlstads Mekaniska Werkstad, Swedish Power _ 
Association, and Hydro - Energy Systems, Inc. (Neyrpic). Mention of these 
_ companies in no way implies endorsement or recommendation by the Water 
and Power Resources Service or the United States Government for these = 
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soivie. EQUATIONS | WITH UNKNOWN 


Analyzing: water distribution networks (WDNs) by computer programs has . 
been carried out for more than a decade [see Shamir for a bibliography updated 
to June 1973 (5)]. The main theme has been to try to help the hydraulic engineer 
in solving those problems which are difficult to solve by hand. However, only 
a few computer programs have been Seeeees to help in the daily operation 
an easy-to-use interactive WDN analyzer seems to called for, 
To popularize the usage | of such analyzers, ‘Programs sh should be developed 
which could be run simply and relatively inexpensively on minicomputers. One 7 
3 computer sysiem developed along these lines in IBM (International Business 
y | Machines) Israel is the Water Distribution Network Analyzer (7) which is based | 
the hydraulic network solver of Brailovsky and Rodek (1) and extensions 
es Eppa and Fowler (3) observed that loop equations tend to consume less computer 
“memory (and probably less ‘computer time) than node equations. Studies made 
in Israel indicate that the number of loops in irrigation water distribution networks 
almost never exceeds 30 (usually it is less than 15) while the number of pipes — 
may be a few hundred, justifying Epp and Fowler’s observation. Shamir and 
-Howard (6) generalized the node equation formulation to deal with different 


types and combinations of unknowns. The unknowns may contain heads, water 
demands, and pipe characteristics. In (6) Shamir and Howard formulate oi 


4 


rules for the proper selection of unknowns. 


_ Inthis paper a method is described for incorporating unknown pipe characteris- _ 
_ tics into loop oriented hydraulic network solvers. A pipe with unknown charac- - 
_ teristics may have either a Positive resistance (indicating that some regulating | 

valve is to be partially closed) or a negative resistance (indicating that a proper 


IBM Israel Scientific Center, The Technion City, Haifa, Israel. 
7IBM Research, 5600 Cottle Rd., San Jose, Calif. 95193. (On sabbatical from IBM 

_ Note.—Discussion open until February 1, 1982. To extend the closing date one month, 

a written request must be filed with the Manager of Technical and Professional Publications, — 
ASCE. Manuscript was submitted for review for possible publication on December 3, 
1980. This paper is part of the Journal of the Hydraulics Division, Proceedings of the _ 
American of Civil Vol. 107, No. HY9, 1981 
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- loop equations with head generators is explained. The intuitive meaning of the 
constituents of the mathematical scheme is considered. Section 4 is devoted — 

to an informal description of an algorithm for automatically constructing and 

solving the set of equations. In Section 5, some of the many possible applications a - 

of head generators are considered. In Section 6 a sufficient condition for the e 

To those of the readers which are not interested in the details of the n attic 


of solving hydraulic problem with head generators we suggest reading Sections P 


on the loop ¢ equations is illustrated. In Section 3, Newion's s method for solving 


Let us consider a network with e e edges (pipes) and n nodes, some of which oa 


2g water sources. (Throughout this paper, only connected networks are consid- 


Mal 


/ eg + 
Yo add FIG. 1.—Example of Network — 


and v, are water sources. Let 2. (i = |, 2, ... €) denote the flow through 
= pipes and let D, (j = 1, 2, ..., m) denote ‘the water demands at the j-th 


node. Then the flow must obey two » laws. 


4 


terms of the Q,’s and D,’s poe n 1 equations (one for each node). Since the — 
_ network is considered to ‘be a closed system, the amount of flow entering into 
= emanating from the network is zero. Therefore, one of the previously mentioned 
n equations eae on the others and thus, there. are n — | independent equations 
_ Conservation of Energy. —Each hydraulic element in the network has its — 
7 on the head loss. Denoting the head loss by Ah the. effect oh! such an element | 


booster must be installed). We shall call these new objects—head generators. 
Davis’ pressure reducing «sda 
— 
Denoting by Ah, the head lo d 


we 
losses around every closed loop in _ network must be equal to zero. If If there 
exists more t than one water source (such as in Fig. d then fictitious | oe 
« connecting on one of them to all the others m may be added to form fictitious loops. — 
In Fig. 1, e, is fictitious. The head loss in a fictitious edge is simply the difference 
in the altitudes of the two endpoint nodes. A network may contain many distinct 
: loops. For example, the network in Fig. 1 contains 3 real (nonfictitious) loops: 
(€4,€5,€ and and three additional loops which 
; contain the fictitious edge e,. Energy should be conserved in each of them. 
_ However, it may be shown that not all of the energy equations are independent. 
In fact, a set of e — n+ 1 loops may be chosen such that the corresponding © 
equations are both independent and imply all the others. (If there are s ‘sources 
s — | additional loops must 
_ The standard technique [e. g.. see (2)] to construct such an independent set = 
of loops is to find a ‘spanning tree (a connected loop-less subnetwork which 
contains all the nodes) and then use the loops defined by the non-tree edges. 
For the network of Fig. 1, a possible spanning tree is the one drawn in Fig. 
_ 2. The node v, is the root (it is the only node into which no edge enters). 7 


that the edges are directed from the root to the leaves. A nontree edge 


rape 


> 


such as e, closes a loop, i.e. set of 


the procedure i is called a fundamental set of loops. the 
-. Let us choose a direction for the nontree edges of the network in Fig. a 


a 
as follows: e, is from v,-v2, e, is directed from and e, is 


+ + + L204) = 0 ae) 
z(v,) and z(v,) = of v, Vas respectively. (Notice that in 
the foregoing equations, only + is used. For explanation see Section 4.) 1 © 
_ To summarize, there are (n — 1) + (e—n+1)+(s— l)=e+(s— 1) equa- 
tions. The unknowns are the flows in the pipes, including the s — | fictitious — 


4 Q’" is a Flow Function which may be Nonzero only, on a ‘Spanning Tree of 


| 
7 
(Notice 
4 
q 
| 


ke 


the Network. ~~, is responsible for conservation of flow and thus it may be 
thought of as carrying the water from a source to the nodes with nonzero — _ 
water demand. Notice that once a spanning tree is chosen, the Q’ flow function _ 
is determined. Thus, for the spanning t! tree in Fig. 2, Q’ is given in the the following — é 


ae of the flow function Q’ for the spanning tree of Fig. “se pvaleotind 


o° is Circular Flow Fenetion.—It may be described ‘ts a set of circular flows, 
one for each loop of the fundamental set of loops. This type of flow is responsible 
for the energy equations. The circular flow variables are denoted by q,, q2, ... 
in which qi is the circular flow in the i-th loop. As explained in the foregoing, 


as will designate the circular flow variables for the loops e>), 


+ Oh +41 +42) + @.)=%- eee 


Ah +4, +q2)+Ahg + + Ah, | 


bh, (Qi + 43) + (QS + +45) + 


computational advantage of decomposing the flow into and is” 
that Q’ is easy to compute and must never be updated once the water demands ~ 
are known, while Q° depends only on e — n+ variables, a number which 
in practice | is much smaller than both e and n. — a sei 
f Till now the unknowns were the flows through the pipes. We wish to ) extend 7 
: ‘the model and allow the pipe characteristics to be unknown. The question may 
be posed as follows: given an edge, e,, what should the additional head loss a 
¥ R, in e, be so that a certain hydraulic effect would be achieved? Notice that 

is a new unknown and that Q, —the flow through e, unknown. 
| Thus, the head loss in e, will be R, + Ah,(Q,). At the present we are interested 
in the ability to obtain a required pressure, P,, in some node v,. In general _ 
there may be several head generators, R,, ..., R,, on r of the edges and the 7 
_ requirements would be to achieve certain pressures in r nodes. In Section 3S 

we show how to use head generators to solve other types of problems. ~e ae 

a It is obvious that the ordinary loop equations are not directly capable 

model head generators. New variables as well as new equations are required. - 
‘Let us consider an omagin first. Assume that a | head generator is put on the 


| 
q 
| 


EQUATIONS 


tes e, of ‘Fig. 1, ead that a pressure, ‘P,, is required at node v,. A new | 
variable R, is used to describe the effect of the head generator. " The bes two 
equations of (2) are not changed. The third equation is replaced by: 
© express the requirement for a fixed pressure in v, an additional equation 
must be introduced: using the spanning tree of Fig. 2 again, the path we be 
the root to v, is and thus we get: 
(QT + 95) + (QF + + 93) + (OF +4, 
+ On Ah. +q2.)+4 + Ah + q2)- ad’ [z(v,) z(v,)) = Qd) 
that the new variable in the now Bq. 24 since 
the path from the root to v, does not aeons the edge on which a head ‘generator 
_ Concerning the solvability of a system of loop equations i in the presence 


of head head generators, three cases may ari arise: cave wi coat 


l. The system of equations has no solution or many solutions. _ 
2 ‘The term R, is positive, meaning that the resistance of the pipe should 


- 


Newton’ s MetHoo For Sotvinc Loop Equations with Heap Generators 


- Solving loop equations is usually done by Newton’ s method (3, 5). 7 To pon 
‘how head generators come into the picture, let us consider first the case where _ 
‘there are no head generators at all. Let Aq be a vector of corrections to the 
circular flow (4q, is a correction to the circular flow q,); J—the Jacobian of © 
‘the system of equations (to be explained later) and h—a column vector whose | 
i th element is th the total head loss in the i ‘i-th loop. T Then, the ‘system to be | 


The content of the matrix J is: a i is the derivative with respect to the 
iow of the head loss in the i-th loop; and (2) J, (i #~ j) is the derivative 
with respect t to the flow of the head loss in the edges common to the i-th 
and j-th loops. Notice that the matrix J is symmetric. Once J and hare available, § 
the system (3) is solved, the Aq are used to update the circular flows ae 
‘new J and h are computed. This is done again and again until the head loss — 4 
Now consider the case when head generators sdo exist. The system of equations Lo 


| to be solved is. is changed into: doide eg 


| 


| 

| 
- 
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 - . The variables are of two types: Aq—as in (3), -" R—a vector of variables, ~ 
iy one per each head generator. 
2. The right hand side also comprises two parts: h' is identical to the h vector ; 
7 of (3) and h* i is the difference between the actual and the required 
a The submatrix A is identical to the Jacobian of a standard network which © 
4, The submatrix B has a column for every loop and a row for every node 
of fixed pressure. For the j-th loop and the i-th node with fixed pressure, = 
; ‘=B, .) is the derivative of the head loss in the edges common to the i -th loop 


and the path in the spanning tree connecting the root to the i-th node with r. 


_ 5. The submatrices C' and c are combinatorial: they contain only zeros 


and ones (the reason for the nonappearance of minus signs is given in Section 


“dl to the loops: an entry of C' is ; equal to | if the j- -th reduction valve is located 
“that an edge which belongs to the corresponding loop. C* is similar to y except 


our using the notation 


Ah; to denote the te devatve “4 
“of 4h,(Q) with respect to the Jacobianiss = 


1 


Ahi Ah! + A 
Notice that the Jacobian is not symmetric. sabe 4 
~ Now we describe an algorithm for constructing e set of ‘equations, seen 


_ Aw Atconitam For Constructing AND Soivinc Loop Eauarions 


a Constructing Depth First Search Tree.—A depth first search tree is a spanning ~ 
tree which is obtained by going as “‘deep’’ as possible: starting at the root, 
try to arrive at a new node; continue from there to a new node; continue 
going deeper and deeper until no new node can be reached; then backtrack — 
and try to emanate again from the node visited before the one exhausted; if | 
necessary, backtrack again. Thus, for the network in Fig. |, the spanning tree 
of Fig. 2 is a depth first sone tree: Start at v,. Traverse e, and arrive . 
Va. . From there, to v, and v, .. From v, there is no way to arrive at a new 
edge since leads Us | to v, which has been visited already» (thus, 


wre) 


| Aa 
| 
pe 
le 
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a nontree edge). Backtrack to v.. Traverse towards Ves then to v,, then 
oy The ede es leads | us back to v, and thus, e, is a nontree edge. Backtrack ¢ 
.. The edge e, is again a nontree edge. Then backtrack to v,, Ys, and : 
v, ant the algorithm terminates. Notice that the algorithm does not try to minimize 
= path lengths. However, by Brailovsky and Rodeh’s construction, there is no 
need to do so. For implementation details see (1). 
- . Besides being very efficient, an important property of depth first search trees _ 
i that a. nontree edge always connects a node with | one of its ancestors in 
; ‘the tree. To show this property, the network of Fig. 1 is redrawn. “Notice 
P the direction chosen for the tree edges is up, while that of the ‘nontree 
edgesis down (thisis always possible) 
> Every nontree edge closes a loop. Traversing such a loop i is done by following | 


aap: edge is met. This explains why “we never have to traverse an n edge 


in a reverse Maection. Even more interesting is the fact that the intersection — 


zhi 


gee 


3.—Network of Fig. 1 
of every two suck loops is a path in a tree connecting a aan ate to one of its 
_offsprings. Thus, traversing the intersection of every two loops is also done 
7 without using an edge in the reverse direction. Now it should become <n 
why only *+’ signs and no * -’ signs appear in Eqs. 1-2. 
generators do not any significant complication. Ifa head 
_ generator is located on an edge e, then all loops passing through it must be 
iene. In our case the head generator is located on e, and thus the only loop 
=. which passes through it is (e, 5355 s€ e,). Again, the choice of edge — 
ensure that only 0 and | (and never i will appear in the submatrices C' 
fay Representing Loops, Paths and Intersections.—To take advantage of the above _ 
_ properties of depth first search trees, the loops, paths and their intersections 


| 
| 


oF 4. The intersection of a loop ae a path— by their eapemaeet common ‘me 
(the lowermost node is the lowermost node of the loop). 
that the intersection of oe paths is not used and thus ‘not be 
4 _ Noticing that the submatrices C' and C® are static (they do not depend in 
the current valves of the variables), the following matrix L may be /Precomputed 
which L" contains the uppermost nodes to the i-th; 
“and j- th loop positions (i, and (j,i), respectively which i > L” contains 
the uppermost nodes commen to loops and paths. 


For out runnin exam le, the matrix L is equal to: 
g p q 


The crucial point is that although the ‘matrix constrection ‘may 
be hard, it is very efficient computationally [see (1)]. = 


Head Around Loops and Paths. the tree 


& 


of s(v) for every node, obtained: 
Wy Bh, + Ah, + Ah, + Ah, 
Bh, + Bh, + bho + Ry 


= 
Ah, + dh, 
Ah, + Ah, + Ah, 


a + Ah, + Ah, + Ah, + Ah, 


c: The foregoing data may be computed by a single ead of each edge. Now 
. computing the head loss along a loop defined by an edge e, connecting v, to 
2 v, is simply s(v,) — 5(¥) + Mh,. | For : a path « connecting the root to a node 
Y v, the sum of the headloss is even easier; it is equal to s(v,). Mere? 
“Computing Derivative of Head Losses Along Loops and Intersections. —This 
os is identical to that of computing head losses along paths, except that 
Bex. the endpoints are found in the matrix L, and that derivatives rather than the = ; 


| 
of the head losses (including the head losses due to the existence of the head x : 
to Fig. 3 the following values 
~~ 


“LOOP EQUA EQUATIONS 

head losses themselves are ve involved (thus, the he head generat eile do not appear). 
— Solving Equations and Updating Circular Flow. —In | Section 4. 2-4. 3a method 
for constructing the system | (4) of « equations thas been p proposed. If the head 
% _ losses along loops are negligible and the head at the nodes with fixed pressure § 
is close to that required then the computation is completed. Otherwise the ] 

corrections Aq must be distributed along the edges of the cornenention circuits. 

_ The nice properties of the depth-first search trees may be used again. For 
; every nontree edge, the only entry of the vector of corrections which —_ 
to . that edge is the one which corresponds to the loop defined by the edge. 

_ Thus, referring to Fig. 3, the correction of e,,e,, and e, is Ag,, Aq,, and 

As for the tree edges, proceed as follows: Assign a value d(v) (initially zero) © 
to every node; add Aq, to ‘the tail and - —Aq, to the head of every nontree 
edge with correction Aq,- For the network in Fig. 3 a values Sencar S are given 
. by the va values of d(v) for every y node, v: ais 


hes 


Now scan ‘the t tree edges in a top op down | manner, scanning an ped “«< if 4 

: the edges emanating from its upper node have been scanned already, and 


accumulate the values of the d(v)’s. The final correction of the circular flow 
in an edge is equal to the final d(v) where v is its uppernode. For the pecetng > 
ex cample we we get the accumulated values of d, for all nodes: on ples a 


Aq 


Again, this method for undetios the flow the edges is computationally 
since it is very efficient. appear in (1). 
aad In this section we consider several practical ste of head generators. 
_ One of the simplest cases is when the given network contains a well and a 
reservoir. The target is to fill the reservoir as fast as possible while still keeping — 
the pressure at the valves sufficiently high. To this effect a - regulating valve 


j 
a 
3 


assume that the altitudes. and water demands at the nodes are as 


meters) ~ 70 0 32 30 35. 

water demand» source reservoir — 150m 2 /h 

Also assume that the regulating valve is on pipe e,, the required 
is 30 m and the pipe parameters are: 

diameter (in millimeters) 150 100 


are obtained: (1) The flow in e, is 220. 31 m n?/h and in it 

_ namely, water is indeed entering into the reservoir; and (2) the regulation valve | 
should create a head loss of 34.1 m. (These results were obtained using the 

- experimental system developed in the IBM Israel Scientific Center. Stet ro 
Now assume that a discharge of 180 m th rather | than 150 m*/his required 

in v,. In such a case the source, ‘¥,, Cannot supply all the demand in the sp 
rear pressure. To obtain a pressure of 30 m at v, a booster must be installed 

on e,. It must yield a head of 47.1 m in a discharge of 58 m m’/h. This time 

: will supply only 122 m 


As a third assume that that a with the following c 
isi installed on é@,: 
efficiency, as a percentage 65 7 5 
a water demand of 180 m */h source v, supplies 99. 1 m ‘the reservoir 
y — 80.9 m */h; the pomp works with 60. 2% ecrey creating a pressure 


_ consumption per cubic meter is increased from 0.226 kW-0. 256 kw but the 
rs total amount of water passing through the pump is lower. To summarize this 
- case, a method to save energy wy? a a regulating valve has been demonstrated. ce 
st 


be ‘inspected with « care. Consider our network again, and assume that pipe e, 
is closed and the water demand at v, is 90 m’/h. Then the pressure at A 7 


Fi 
foll 
g 
| 
: 
™ is introduced on e, to lower the pressure at v, to 30 m the results are: ve 
supplies 122 m’*/h; v, supplies 58 m*/h; the pump works with 67.6% efficiency 
a3 creating a pressure of 62.7 m consuming only 14.9 kW. The regulating valve — - 
q 


unexpected price: a booster, creating a head of 30.6 m for a discharge of Te 
-m’/h from v,—v, must be installed. In this way the flow through e, is increased, 


the head losses increase also and the pressure at v, is decreased. We doubt 


for the “efficiency of the proposed algorithm, all the ‘solations quoted. 
_ above require at most 5 iterations (some took only 4) to achieve a maximal ; 
head loss imbalance around loops of less than 0.01 m. tsod ban 


On Sorvasiuty of System oF Equations In Presence OF Heap 7 


It is well that in the of head generators the of equations. 


“the corresponding set of Eqs. 4 may have no solution or many solutions. The - 
network of Fig. 4, in which v, and v, are nodes with fixed pressure, v, 4 ; 

_ the water source and e, and e, are ofess with head generators, is in general +, 
not solvable: ‘the two head | generators” have an effect of one and thus, there mi 


and nodes with fixed pressure in the network, a sufficient condition is developed. — a 
_ First, notice that we may assume that the network contains only one water 
source. Otherwise choose one of the ‘several water sources: as a a reference node, 


__ pressure. . An obvious ‘necessary condition for solvability of (4) for i Peon 
. pipe characteristics is that the number of head generators equals the number ¢ 
_ of nodes with fixed pressure | (not counting the reference node, , although | the 
a _ A sufficient condition for the existence of pipe characteristics for which (4) 


has a unique solution is the existence of a connected subnetwork G’ Go of ‘the 
_ given network G such that: 


wh All head generators, nodes with fi fixed pressure an and the d the unique \ water source ec 
belong to G’. 


. For every two nodes with fixed pressure (or a node with fixed pressure — 


and the water source) the unique path in G’  E them contains at least 


on! { iss 


i 
| 
pressure requirements). This type of unsolvability of (4) is due to the topological [i 
«FIG. 4.—Two Head Generators with Effect of One 
| Me Dy Cases Tesistalice Of Wich 
a 


‘The proof that these conditions indeed sufficient is given in 
_ Notice that the network of Fig. 4 does not satisfy the aforementioned sufficient 4 
c conditions. If the water source were v, rather than v, then the required pressures" SS 


 Anetwork which topologically seems all right may be problematic with respect — 
‘to certain data. Consider, e.g., the network of Fig. 5. Assume that v, is a 
_ water source, v, is a node with a required pressure P, all the edges are identical 
and e, contains a head generator. If the problem has a solution with a nonzero 
_ flow through the pipe e,, then the reverse flow is an alternative solution, as 
can easily be seen by symmetry. Therefore, two solutions are possible. bP +O 
If the altitudes of all the nodes are zero and - > 0 then the problem has 
om no solution of all. The reason is that the head generator must have a nonzero 


effect. Without loss of generality we may assume that the flow ‘through o 
goes from v,-v,. The direction of flow must be as in Fig. 6. Assuming an 


_ water demand D D,zaty and denoting the flow through pipes b Bs we — 
thus + - - D,)= )+ Ah(Q,- D,). 


Since 4h is a monotonic function this ‘implies that Q,= = 0.3 therefore Q. 


oni extensive with head | generators has been conducted 
in Israel during 1978-80. It was found that the method is very useful for both 
- designing and operating water distribution networks. The computer time needed 
to construct and solve the system of equations is negligible. In fact, even large — 
; networks with about 300 pipes and 25 loops, water sources and head generators 
can be solved efficiently on a portable computer such as the IBM/5110. ne “ 
is believed that this would have been impossible if node equations were involved, 
due to storage limitations. An easy-to-use interface was inevitable if these ideas 


at v, cannot t be 


and indeed, many hydraulic engineers, students and even farm managers have 


=f 
> q 
OT 

4 
| 
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We wish to acknowledge fruitful conversations with hea great 
help i in carrying out experiments given to us by J. Heller of the Israel Ministry 

; of Agriculture, and Z. Shavit of the Israel Water Workers Association. ty 


1.—Proor of Surricient Conpition OF Section 6 wr 
 «G'isa loopless subgraph of th the given network, G, | all head 
_ generators, nodes with fixed pressure and the reference node. G’ may be 
augmented by nodes and edges to form a spanning tree, T, rooted at the reference 
“node. The head generators may be renumbered eneneting to their distance from 3 ; 
= root such that R, is the closest to the root, R, is next to closest, etc. 
7 If two head generators have equal distances men the root, the order ‘between 
a. every node with fixed pressure we associate the head generator closest 
to it on the way to the root. Condition (3) implies that this association is biunique. 
Otherwise two nodes with fixed pressure would be associated with the same 
head generator, which implies that the path connecting them in the tree does 
- not contain any head generator. Let us renumber the nodes with fixed p 


- such that the i-th node (denoted by u, ,)is associated with the i-th head generator | 


Cs 


_ Referring | to Eqs. 4 | constructed for the spanning tree, T, the matrix C*h 


Thus, the matrix c ~ iS nonsingular. It ‘has an additional interesting property: 
if a) #0 then the first j entries of the i- th row are all equal to the first 
of entries of the j- -th row. Thus, if k is the maximal column number such that 
= # 0 and k < i, then by subtracting the k-th row from the i-th row, wall 
i-th row will contain a nonzero element only in positionC;. 
consider some initial flow which is zero in all the ‘nontree “edges. 


' - Perturb Q” by adding a circular flow of 1 m >/h around each of ‘the loops, 
Y to obtain a flow Q. ‘Choose pipe parameters such that the resistance of each 


nontree edge is very high while those of all the tree edges are very low. The 
matrix A of (4) will have big elements along the diagonal and small elements — 
_ elsewhere. The elements of the matrix B will also be small. Thus, by phenere ll 4 
rows to transform C* into a unit matrix, the elements of B remain small. Then 
_ the rows of C* may be used to transform C' into a matrix of zeros only, 
leaving the diagonal of A large. This shows that the Jacobian of of (4) sion Ard 
made nonsingular, and therefore, that (4) is solvable. 
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submatrices of of of equations w when head 
generators are present; 


an auxiliary vector to compute to flow 3 ls 


number of edges in network; 
of head losses around loops; 


Pp 
Jacobian of standard problem (with no iis generators); 
~ matrix representing intersections of loops and paths; _ 
submatricesofL; 


= number of nodes in network; 


ot 
vector of circular flow; 
vector of flows which is nonzero only o} on spanning tree; 

number of sources; 


- sum of head | ~ on path from root tov 
Ah; = de erivative of head loss in e,; 


vector of head losses; and 


vee 
_ vector of corrections to circ to circular flow. dhe 


7 
zz. 
4q 
= _ 
- 
| 
uf | = 


‘ng the quant itics; €, in the 
8 Srupy oF F OF FLoop Wave PROPAGATION 


‘The: one- dimensional, ‘varied, unsteady flow eq equations, commonly 

_ referred to as Saint-Venant equations, form the basis of study of flood wave 

. motion in open channels. Recent studies in this area were generally concerned — 
with finding new methods for solution of these equations. A number of numerical — 
schemes have been proposed based on finite difference methods. These schemes — 

can be broadly | classified | into the characteristic 1 methods i ‘in which the integration 


basic equations are directly expressed in finite difference form (1,5,7,13). The 
direct methods can again be classified into explicit (5,7,10,12) and implicit finite 
difference methods (1,7,10,13). The explicit methods which advance the solution | 

, by simple explicit algebraic expressions suffer from the disadvantage of restricted — 
stability. In the | implicit method, ‘unconditional stability is obtained at the cost 

4 Most of the earlier studies on numerical solution of flood wave problem 
q either deal with the task of developing new numerical schemes, improving old — 
Remeccite or applying these methods to specific field problems. There have been 


- an initial uniform flow depth of 0.3 ft (0.0915 m). A ‘ee wave of =i 
sec duration was taken as the inflow hydrograph. Using the characteristic method 4 
7 on a rectangular grid, they studied the effect of wave amplitude, Manning’ . 
_ nand bed slope on subsidence of the wave. While these studies g give an indication 6 
7 of the nature af damping of a flood wave in a prismatic channel, the studies — 


; 7 0.002 ft (0.610 mm) to 0.02 ft (6.1 mm) forms a very small fraction of the 
— initial flow depth of 0.3 ft (0.0915 m). A qualitative study of subsidence in 
- prismatic channels for such small amplitude waves was also made by Henderson © 


"Asst. Prof., Civ. Engrg. Dept., Indian Inst. of Science, Bangalore, India 560012. 
Research Scholar, Civ. Engrg. ‘Dept., Indian Inst. of Science, Bangalore, India 560012. 
_ Note.—Discussion open until February 1, 1982. To extend the closing date one month, | 
a written request must be filed with the Manager of Technical and Professional Publications, — ; 
ASCE. Manuscript was submitted for review for possible publication on December 4, 
1980. This paper is part of the Journal of the Hydraulics Division, Proceedings of 
- American Society of Civil Engineers, ©ASCE, Vol. 107, No. HY9, September, 1981. 
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In order to obtain a a generalized picture of the « damping 2 of a a flood wave 
_in a prismatic channel, it would be better to convert the problem into a 
nondimensional form and make a parametric study to investigate the effects “4 
of flow parameters, wave parameters and channel ‘shape parameters. The ese = 


The present paper presents a parametric study of flood wave propagation ~ 
problem based on numerical solution of the nondimensionalized unsteady fl flow 
equations. ‘The ¢ propagation of a sinusoidal flood wave in a prismatic channel 
_ is studied for uniform initial flow. The governing parameters, namely the initial 
uniform flow Froude number, the aig parameters and the channel shape 


parameters are varied over a wide range. 


= The one dimensional equation | of continuity and motion for unsteady chann = 


. 
in n which x along = time; y = dept th of flow; 


_ ‘velocity of flow: a = cross-sectional area of flow; w = top width of flow; q 
So = bed ays, S, = friction er and g = gravitational acceleration. if 


4 


3 


in which r = hydraulic radius and subscript o refers to initial uniform flow 
conditions. Eqs. | and 2 are nondimensionalized to enable a generalized parametric — 
study, as follows. The depth y is nondimensionalized with respect to y ,, initial : 
uniform flow depth; all cross-sectional dimensions are nondimensionalized with 
_anifor to y,; the longitudinal distance along the channel x is nondimensionalized _ 
— with respect to/, = y,/S,; the velocity v is nondimensionalized with respect 
to v,, initial uniform flow velocity; time ¢ is nondimensionalized with — 
= | With nondimensionalization, Eqs. and 2 become, 


| 
| 
- _ study for unsteady flows was made by Sakkas and Strelkoff (11) for the dam 
q 
ay av av | 
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FLOOD WAVE PR PROPAGATION 


depth corresponding | to a depth y, given by a 


is the uniform n flow Froude number given by by abd bas 
, and are of Y given by za 
ea barton ‘ed banies du 
‘For re and trapezoidal channel shapes considered in this study, d 
, and , are functions of B, and Z in which B, = y,/b is the inverse of 
a nondimensional bed width with b = channel bed width and Z is the side 
slope of the channel (Z = Ofor rectangular channels). 
= The p present studies are concerned with the propagation of a a sinusoidal flood 
wa wave along a prismatic channel in which the initial flow is uniform. At the 
deft boundary, X = 0, a sinusoidal inflow hydrograph starting at T = 0 is 
considered. The wave introduced atx=Oisoftheform 
wale 
; O<t<t,y=y,; t 
in which Vw = wave amplitude; and ¢,, = duration over which the wave wend 
= 0. ‘Nondimensionalization of Eq. 10 leads to 
< 


COS 


<T,Y=1; T 
a 


4.4 


ts 


1; and V = 1 may re treated a as 


hei othe sk solution region. lo bos # woled O = 


pare parameters the problem defined by Eqs. 4, are the 
- initia flow Froude number F., nondimensional wave amplitude Y,,, nondimen- 


sional wave duration 7,, and the channel shape parameters. For rectangular 
and trapezoidal ct channels, the > shape parameters a are B, and Z. These parameters 


: are varied over a . wide range systematically and their effects on subsidence 
I the wave is studied. The initial uniform flow number, F, is varied over 


o 


_ 7 _ in which capital letters denote the nondimensional quantities with the correspond- ; 
| 
¢ 
7 
| 
| 


range 

0.5-3.0; B, is varied from 0-1 to wide | channels to 
very narrow channels; Z is varied from 0-3. As the study is concerned with % 
subsidence of flood waves, the wave parameter values for different cases are A 

chosen such that shock formation does not occur, 
. An irregular characteristic grid method has been used to get the solution. 

- and Woolhiser (7) has compared this method with other numerical methods 
and claim this method to be the most accurate one. A disadvantage of this 
w_— is the need to interpolate the results for depth and velocity at fixed © 
sections or fixed time as the solution is obtained at grid points governed by 

the characteristics. In the present computations, interpolation was done contin- — 


_ uously as the computations progress. Two cases studied by Mozayeny and Song 
Owe were 'e solved using the characteristic grid method after n ‘nondimensionalization 


obtained by “+ emet and Song (9) using characteristic method as a rectangular _ 

«grid compared well. The subsidence indicated by Mozayeny and Song 

- larger with respect to the results from the characteristic grid method. | 

Computations were on the DEC 1090 digital computer for cases 

f= in all and for each case, computations were made for over 50,000 grid points — 

j in the XT plane. Effects of F,, Y,, T,, B,, and Z on the modification of — 

_ Stage hydrographs, discharge hydrographs, instantaneous wave profiles, stage 

ad discharge subsidence rates and the speed of the wave peak were studied 
in detail. The results are presented in a summarized form here > laying stress” 


Stage and Discharge Hydrographs.—Typical results of the modification of stage 
and discharge hydrographs are presented in Figs. 1 and 2 respectively, for a 4 
wide channel (B, = 0) with F, = 0.3, Y, = 0.5 and ae =3. .0. The nondimensional — 


— is defined by 


(13) 


BAS 


ie. 


in which discharge at at any location and and q 


of the wave with distance, which was observed in varying degrees for all the 
cases. The rate of rise of the hydrograph is found to be greater than the rate 
_ of recession as typical of a flood wave. The tail end of the discharge hydrograph 
xX =0 dips below 0, a trend observed by Mozayeny and Song (9) also. 
¥ This is because in the recession stage of a flood, the same stage yields lower g 
_ discharge due to the nature of water surface profile along the channel. irene 
Effect of Froude Number on Subsidence. —Subsidence of mo is studied in if 
gives the local wave amplitude normalized 


= 


q 
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FIG. 1.—Modification of Stage Hydrograph with Distance, | 
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arge,| 


_ with respect to the initial wave amplitude at X = 0. Subsitenee of ss. 
is considered in terms of an analogous quantity, the rise 


PEt 

Og 


Q* with X for different F, values are presented a 


= 


t 


~ 


and F, > 0.8, local supercritical flow might occur during the time of peak : 
of the wave. Thus, computations were restricted to F, < 0.8. The significant _— 


| 
tz: 
4.—Subsidence of Relative Discharge Rise—Effect of Froude Number 
: 2.055. 2.057. and 2.058. respectively. Computations showed that for Y. = 0.5 
| 
| 


FLOOD | WAVE PROPAGATION 


damping of the wave in “prismatic channels | is clearly out by Figs. 
and 4. The wave amplitude at X = ‘i and 10 are 59% and A% respectively — 


7 being more at lower F , values. This 1 may be explained as follows. If we consider : 
a specified initial uniform flow depth y , and chennal bed slope S » from Manning's 
formula an and the definition | of » we have 


The length scale 1, is constant as y, and S, are kept constant. The time scale 


_ Thus, under the conditions specified, an increase of F, denotes a poor 
of channel roughness and time scale. The decrease of time scale implies a 
flood of shorter duration for a constant nondimensional wave duration, T 24 
_ Computations indicate that a a a reduction | in the wave duration would result in 
a larger subsidence of the wave. On the other hand, the decrease of n associated q 
: with the increase of F, would indicate a lesser subsidence of the wave. The 

- fact that the subsidence i is lesser for higher F , values indicate the greater influence 
_ of the second factor, namely the decrease i in channel roughness. Computations | 
for different combinations of Y , and consistently confirm this observation. 
Zz Mozayeny and Song (9) reported an exponential decay for the wave amplitude _ 

. 4 with distance. The present results show that subsidence in the initial reaches 


. 


| 
| 
: 
var 


is in an exponential re relationship such as 


i 

as the \ form modified with distance, a ‘constant “exponent 
_ does not describe the subsidence rate for the entire reach. The results show 
that corresponding to a greater subsidence and therefore greater change in the 


wave form, the limiting X up to which K Temains constant is less for lower 


| 


‘*F, values. For the results given in Fig. 3, K value for the initial reach is found 
to vary with F., from 0.104 up to X = 5 for F, = 0.1 to 0.075 up to X 
= 8 for F, = 0.7. The channel reach considered by Mozayeny and Song 
to X less than 0.35 and thus a constant was obtained for the 


po 


Di t x 


6—Si of Relative Wave Amplitude—Effect of Wave jo 
A plot of Y*® against F, with X as the third parameter indicated that Y* 


varied practically linearly with F. for a constant X. Thus, an equation was” 


0. 5 and a a’ 0. Eq. 19 can form the basis for a generalized prediction of _— 
subsidence in wide channels where the coefficients, C, and C,, are a | 

Effect of Wave Amplitude on Subsidence. —The Variations of and Qt with 
X for different Y, values for F, = 0.3 and T, = 3.0 are presented in Figs. 

6 and 7 respectively, for wide channels. The value of Q.,,,x (0) for Y, = 0.1, 
0.5, 1.0, and 1.5 are 1.186, 2.055, 3. 354, and 5.05, respectively. It is ‘clearly 


| 
— 
| 
fe 
; 
| seen Irom Figs. 0 and / that wave amplitude has some eifect on the relative 
subsidence rate confirming the nonlinearity of the phenomenon. Subsidence 


is relatively lesser for wave 2 ‘this effect is. 
clearly lesser in results of relative discharge 1 rise. In fact beyond a particular 
X value, the Q* values are practically identical for different Y,, values. The 
‘wave amplitude at X = 10 is 48% and 42% of the initial wave amplitude at 
‘Xx = 0 for Y,, = 1.5 and 0.1, respectively. The corresponding discharge rise 
at X = 10 is 39% of the adie discharge rise at X = 0 for both Y, 
Effect of Wave on Subsidence.—The variations of Y* and Q* with 
X for different T,, values for F, = 0.3 and Y, = 0.5 are presented in Figs. — 
8 and 9 respectively, for wide channels. The values of Q,,,, (0) for T,, = 0.5, 7 
1.0, 2.0, and 3.0 are 2.486, 2.274, 2.116, and 2.055, respectively. It is seen 
that wave period has a pronounced effect on the subsidence rate. The wave 
amplitude subsides to 13% of the original amplitude at ie + a for T, = 0.5— 


e rise,Q 
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> 
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while the corresponding ng value i is 60% for T, = =3.0 (Fig. 8). The 1 rate of aiabiiains 
- 4 is very large in the initial reaches for low T,, values and the subsidence rate 
- comes down only after the base of the hydrograph has spread significantly - 
at a sufficiently downstream location increasing the local wave duration. The 
results in Figs. 6-9 indicate that as the bulk of the wave form reduces, either 
through reduction of or the subsidence rate But a variation 
in T,, has a much larger influence than a variation in Y,. 
4 Effect of Shape Parameters on Subsidence.—The variations sof Y* and Qt 
with X in rectangular channels (Z = 0) for different values of B, are shown ; 
in Figs. 10 and 11 for F, = 0.3, ¥, = 0.5, and T,, = 3.0. The values of > 
for B, =0, 0.1, 0.5, ‘and 1. Oare 2.055, 1.969, 826, and 1.767, respectively. 


= results show that there is a Significant effect of the width parameter =. 


| 
aS 


and the length scales are taken as constants, nts, and S, are ire constants. 
_ Thus, for a constant F,, Manning’s m decreases as B, increases. Therefore, a 
when computations are made for a constant F, value, a change in B, also” 
implies a mpg ol in n value as governed by Eq. 20. Thus the results for narrower 
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AIG. 8. of Relative Wave Amplitude—Effect ot Duration 


= «FIG. 9. 9.—Subsidence | of Relative Discharge Rise—Effect of Wave Duration 
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to with, a n value than for 4 


y wide channels. In spite of this, a greater subsidence rate is observed for narrower 

_ channels (Figs. 10 and 11) indicating a very significant effect of the width parameter a 
on the rate of subsidence. The side walls have a pronounced effect on the oa 
subsidence in narrow channels. == 
“| The effect of channel side ¢ slope Zon ae is shown in Fig. 12(a)-12(c) . 


in which Y* variation is given for Z = mm 0, I, a. and 3 for constant B, values. 
‘Wave 
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0.—Subsidence of Relative Wave of Width Parameter 
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FIG. 11 Subsidence of Relative Discharge Rise—Effect of Width Parameter 
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a pm Z as subsidence is pronounced in narrow channels (Fig. 12(c)], while, 
it is insignificant in the case of relatively wide channels [Fig. 12(a)]. For B, 
= 1, the rectangular channel has clearly a greater subsidence while there is 
again only a small influence of Z in the range 1 = Z = en 
For channels, the equation corresponding to Eq. 20i is 


= we consider y,, S,, F,, and Y,, as fixed, the length scale i is constant t while 
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plitude, 
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Relative wave amplitude 
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FIG fa, of Relative Wave of Side Slope 
the ena s n and time scale vary with B, and Z as given —" 21 and a 
_ 22. For a constant B,, n and ¢, both increase with j increase of Z. While increase ; 
- of n will lead to greater subsidence, the increase of ¢, and thus ¢,, (for constant 
T,,) will lead to lesser subsidence. The results in Figs. 12(a)-12(c) show that 
; except in narrow channels, these two factors largely omen each other leading — 
an insignificant influence of Z on subsidence. 
Speed of Travel of Wave Peak.—In studying the propagation of a flood wave, —, 
two principal questions: arise, namely, the magnitude of the wave ve peak at a 


| 
| 
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= particular location and ‘the 1 time at which the peak a: arrives at the location. The 
path of the wave peak in the XT plane was studied for all the computational | 
eases and only a summary of the results are indicated here Computational 
results show that the nondimensional ‘speed of of the wave peak does: ‘not vary 
significantly with F , 
increases as the velocity scale i is directly proportional to wi This greater speed 
of the wave peak corresponds to the lesser subsidence at higher F > Values. 
‘Initial wave amplitude Y,, has a influence on the speed 


significantly a: as seen . earlier, it has only a "small influence on the speed of oe 
wave peak, with the speed being slightly larger for higher 7, values corresponding — 

- to lesser subsidence. The nondimensional speed of the wave peak is found 

. | to be practically independent of the side slope Z except in very narrow channels 
in in which the speed is slower in rectangular channels Corseapending to larger 


a ‘Kinematic wave 


pas 


13.—Path of Wave in XT Plane, F, = 0.5 
‘wave subsidence. The speed of the peak is lesser in in narrower channels. (larger 
corresponding to greater subsidence, 
The computational results were used to compare the speed of the wave peak 


with that of a corresponding | aici wave. The speed (v,,) of a kinematic 


_ Nondimensionalizing Eq. (23 (with Tespect to y o) for wide channels and using 


nen 
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— ndimensional kinematic wave speed can be calculated from Eq. 24 based | 


the speed of the flood wave peak v very well. The conc concept - of a kinematic wave 
- implies that the bed slope term dominates over the other slope terms in the = 
_ equation of motion (4). Therefore we might expect the flood wave to move 
more closely like a kinematic wave at higher F, values and the computation 

ig The path in the XT plane along which the depth is 1% _ above the: ‘initial 


7 characteristic as the 1% increase » Y denotes a larger percentage of the mee 


= An extensive parametric aie on the propagation of a sinusoidal wave in 


open channels leads to the following conclusions: 


There can be significant subsidence in Prismatic open channels. ‘This i is 


Initial ‘uniform flow Froude number F. has some jnflwence on 
with greater subsidence for lower F, values. The relative wave amplitude Y* 
varies linearly with F,. The nondimensional speed of the wave peak does not | 
values” giving slightly higher “subsidence. “This nonlinearity | is “found t o be 
lesser i inthe discharge results. 


4, Wave duration T,, has a pronounced 1 on subsidence with subsidence 
being significantly more for lower T,, values. 


‘significant ‘except in ) narrow channels. The nondimensional speed of the wave | 
peak is also found to be practically independent of Z except in narrow channels. | 
The width parameter, B,, has a significant effect on subsidence, which is least 

in wide channels corresponding to B, = 0. The speed of the wave peak is 
lesser i in narrow channels corresponding to greater subsidence 

ce ‘All computational results show that the ‘relative wave amplitude, y: 

the relative discharge rise, Os. decrease ¢ exponentially with distance in nthe 
initial reaches. This exponential fit will hold good up to a larger distance in _ 
cases where subsidence in lesser, 
7. There is a very good agreement between the speed of the Sa peak 7 
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= area of flow; 


inverse of nondimensional bed width, y 


= initial uniform flow Froude number, Ne d, 7 
acceleration due to gravity; 
exponent in Eq. 18; 


nondimensional discharge, g/q,; 
relative rise ‘tise (Eq. 


@ 
| 
&§ 
— 
| 
aper: 
C= coefficient in Eq. 19; 
| 


qa = to initial flow; 


“ rs = nondimensional t 
T. 


= time scale, / 
= 
WwW = nondimensional top width, w 


* 


x= = distance alongchannel; 
Y = nondimensional depth of flow,y/y,; 

Y, = initial wave amplitude, y,/y,; 


Z side slope of trapezoidal channel; 
= function given by Eq. 8; sand 

function given by Eq. 9. 
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DemMaA ND FORECASTING IN ‘WATER 


The ray rapid increase in energy costs in n recent years has focused attention more 
sharply on the potential for energy savings in the operation of water distribution 
networks. The largest single cost element is often the cost of operating the | 
pumping stations to meet consumer demands, in terms of quantity and pressure, 
distributed d throughout the network. A fundamental requirement for Scheduling 7 
pumps in an optimum “way is a forecast of what the consumer demands are 
_ likely to be, say 24 h ahead in time. Provided this forecast is reasonably accurate, 
the reservoir levels can be allowed to rise or fall according to whether the 


_ networks, however, there is is s considerable potential for e energy | saving by petrwn 
- pumping for electrically driven pumping stations with variable speed drives. 
- Such situations require an accurate forecast of the demand variation throughout 
the day. This more difficult problem is examined here. 
This paper is concerned with formulating the h hour by hour domend forecasting 
- problem within the framework of an on-line | pump scheduling ng scheme, and w with 


c 


comparing two well-known statistical forecasting methods to evaluate their 
effectiveness in terms of accuracy and computational feasibility. The computa- - 
a tional results are based on data from the East Worcestershire Waterworks y 
7 company in the United Kingdom (5,6,16). The first method uses Kalman filtering 
techniques with a temperature or rainfall variable controlling the propagation — 4 
: of the residual component 0 of demand. The second method is based on the | 
Karhunen-Loéve spectral expansion and has the advantage that it uses only = 
past load data and does not require any meteorological inputs. © Arts Big hr 
Load forecasting is still considered more of an art than a science (9, 14). 
For this reason many distribution industries, such as electricity, gas, and water 
supply, maintain a number of load forecasters who base their prediction on. 
‘experience and insight. Basically they | maintain records of past load behavior, 
'Lect., Civ. Engrg. Dept., Imperial Coll., , London SW7 2AZ, England. 
Note. Discussion open until February 1, 1982. To extend the closing date one month, 
‘a written request must be filed with the Manager of Technical and Professional Publications, 
ASCE. Manuscript was submitted for review for possible publication on June 3, 1980. 
BS This paper is part of the Journal of the Hydraulics Division, Proceedings of the American — 
Society of Civil Engineers, ©ASCE, ‘Vol. 107, No. Septeasber, 1981. ISSN 
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actual demand is below or above forecast, enabling pump energy savings to 
7 " be made by operating the pumps at constant output. For such purposes, the = 
main requirem is an accurate 2 recact man many 
i 
. 
: 


‘together with pas past weather ‘conditions, industrial series, and other factors which — 
are known to significantly affect demand. Based on this information a load 7 
4 prediction is made by comparing the conditions with those of a similar day 
some time in the past, whose characteristics they have recorded on file. The 
resulting forecast is frequently quite adequate. aw aeratl 
The disadvantages with this approach lie in the fact that large amounts of 
data have to be stored, only one or two forecasts can be made daily, and | 
the method is highly subjective and therefore liable to human errors and biases. — 
“Abernatively, this intuitive approach is sometimes advantageous if the complexity 
of the problem is such that human insight may be the best way o! of arriving — 
at a forecast when conditions deviate Significantly from the norm. 
Over the “past few years « a serious attempt has been made to develop 
_ mathematical models which can be implemented in a computer for automatic — 
load forecasting. These methods, which have been well surveyed by Matthewman 


Z and Nicholson (13) and Gupta (5) are of two main types. One method is to 


establish a ‘mathematical. model based on the correlation between load and weather, 


Bo. 4 raised since they require weather data monitoring, which is sometimes 
: ‘not readily available for on-line implementation in the United Kingdom, although 
it is in California. Techniques using only past load data have the advantage © 
: of being much simpler to implement in the sense that data is more readily 


available. Also it is argued that some weather e effects on the load can edly 

_ relatively long time constants, of the order of 24 h, and are therefore easily — 

identifiable and predictable in terms of the most recent past load data. On 

_ the other hand, in cases where weather effects have shorter response times, 

s the order of 2 h-3 h, it seems quite sensible not to reject evident correlation - 

between them and the load pattern. The use of past and present we ather \ would 

-;- to be ve very significant in arriving at better load forecasts ‘for lead times 

of 24 h or more, the argument being that past load data alone does not contain 
; Bed about the present and most recent weather effects. In addition, — 
; _if weather variables are incorporated into the model, then the effects of sudden 
changes in weather conditions can be taken into account and perhaps contingency © 
predictions obtained for different possible weather conditions. 
_ Possibly the widest use of short term load estimation has been in the eet ‘ 
ae industry. Work on load models based on weather correlations has been 
deme by Dryar (4), Davies (2) and Heinemann (10). The basic approach is to 
find a deterministic relationship between the peak load and the average daily _ 
weather effects } considered | significant. Load variations with time of of day are 

"not considered and there is no measure of confidence or uncertainty y associated — 

_ The most important load correlation technique is the method of Farmer and — 
Potton (7) who introduce a probabilistic structure into the model. Essentially, — 

— load observations overa period of about 6 weeks in the past 2 are used t to evaluate 

to model the load using a finite Karhunen-Loeve ‘expansion The parameters 

= the expansion can then be recursively updated to fit the most recent 
| 


observations. This method, which has been used extensively in practical applica- __ 
" tions, has the advantage that it provides a continuous forecast and error standard — 


| 
q 
| 
q | 


the first to suggest a a state space both the effects 
- of weather, temperature, and humidity, as well as the effects of the most recent 
demand observation. The method uses Kalman filtering techniques for providing 
continuous forecast and error standard deviation. Galiana (8) used a similar 
model and extensively tested the method on real data. a6 


In the case of water ‘supply, De Moyer, Gilman, and Goodman Q) performed 

some preliminary work on demand estimation using a Fourier series method 
: _ with an associated accuracy determination. So that no more terms than necessary 

-= calculated for the series, a statistical test was used to indicate the number 

of harmonics which significantly add to the series representation of actual demand — 
“# ‘data. . For summer consumption, maximum daily temperature is used —— 

; ‘ther prediction accuracy. Shimron and Shamir (15) use time series analysis methods 
for forecasting water demands several hours in advance using recorded hourly 
readings. Two basic models are used—one is a ‘“‘means’’ model, which forecasts | 

future demands as weighted averages of selected past data, and the other i. 

“‘additive’’ model, which decomposes the historical data into several basic 


time series, extrapolates each time series ‘separately into the future, and then 
_ adds them together to produce the forecast demands. The disadvantage goood - 
_ both of these methods is that the accuracies obtainable suggest that they are 
inappropriate for use in on-line optimal controhh 
The two primary aims for a good prediction method are accuracy of forecast 
_ and ease of implementation. Exact numerical accuracy de depends to a large extent 7 
on the particular network and data being analyzed, and will also be closely » 
related to the length of forecast required. Shorter prediction times result in 
— accuracy, but if forward predictions in excess of 24 h are required, then 
errors will rapidly become unacceptable unless more sophisticated wangee dl 
are used. On the other hand, emphasis must be placed on n techniques which 


are : readily implementable ot on- -line and will ‘not burden the ; available computer 


excessively. Methods which collect and process data on-line to update the current 
forecast are clearly desirable, as are methods which give a continuous record — 
os of the model’s confidence in the forecast or estimate of the prediction error. 1% 
The two methods outlined in this paper have been specifically designed t to satisfy 
most tof t ‘these tin 


_ Forecasting plays a central role in the on-line scheduling of pumps in water — 
_ Supply networks. The main reason for this is that the savings to be gained 
by hourly “fine tuning’? of pump schedules are limited the with 
_ To fully appreciate the necessity for accurate on-line forecasts, it is necessary _ 
to look more closely at the relationship between the on-line forecasting problem 
and the problem of on-line pump control. They are both computationally | 
_ demanding tasks which should be combined with data acquisition and monitoring 


a water network to to be updated and schedules to to 
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_ deviation. A similar approach is presented by Christiaanse (1) who models hourl . 
q 


q 
q 
| 


cumin system, which measures network data z and c carries out pump on nea ty 
and forecasting calculations, and an actual network consisting of pipes, pumps, 
_ reservoirs, valves, etc. The basis of this interaction is that the computer must © 
construct demand estimates based on the latest observed network parameters, 
— as pressures and flows, since on most water supply networks there is 
no direct on-line metering of consumer demands. When these forecasts have “lt 
‘been made, an optimization calculation is carried out to determine the optimum — 
pump schedule which is then implemented on the network. The important feature 
of Fig. 1 is that the optimal control calculation is based on a forecast which | 7 
in turn copends on actual on-line pressure and flow measurements from the . 
- network. In other words, , the mechanism for feedback control from the network _ 7, 


> 


Consumer 


Water Supply 


ear, 


by the computer to provide a basis for deterministic optimal control of pumping 4 

z Among the state variables on the network, as will be shown in the next section, 
are the residual components of the water demand. Thusthe pumpcontrol problem © 
can be viewed as a state feedback problem in optimal control where there 4 

is uncertainty in both the state propagation and measurement equations. ‘The: g 
problem of estimating the state under such circumstances can be solved using fi 


the Kalman filter approach, which is the first of the techniques considered 
for demand forecasting, and is clearly suitable to on-line implementation due 
to the advances made in the computational aspects of filtering since the 1960s. a 
An alternative to the feedback filtering approach | is the ae expansion 
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two methods of forecasting considered. 
= _ This paper examines both of the forecasting methods to determine which _ 


NETWORKS 


is most suitable for use in an on- on-line c ‘computerized pump scheduling. ‘scheme. _ 
_ Several such schemes are currently in operation, both in the United States | 
_ and Europe, and in most cases forecasting is done on an ad-hoc basis or we 
_ some very simple modelling relationships. Such rudimentary forecasts are not — 
Ea appropriate t to » sophisticated computerized systems which rely on an accurate a 
forecast to generate savings in pumping costs. 
7 5 ‘The water demand at any location on a water distribution network depends — 
on the pattern of water consumption of the individual industrial, commercial, — 
- public, agricultural, and domestic users in a given area. On aggregate, these 
can be assumed to exhibit a cyclical trend which can be modelled as a time 
5 series with inherent uncertainty. Meteorological factors may affect demand in © 
a systematic way, e.g., domestic consumption during the summer can “ie 
affected by high temperatures. is often difficult, however, to to 


period can 1 be modelled as a Fourier Series 
q 


walk 

y,(k) = ar + (a? sini wk + b? cos iw ok) 

‘in which n p> a’ and b? are to be determined; and w, = the natural frequency. . 

_ The residual component can be written as an auto-regressive- moving- average 


(ARMA) model of the form 
ara 


in which | u=a parameter; w = white noise; and m, 
, d, are constants to be determined. The w noise has first and second order 
statistical properties 


expansion of eigenfunctions of the demand data covariance 
Thus for a particular day period k of N discrete intervals, the total water demand 
is represented by the mth sample function q,,(k) of an ensemble of Mdays 
(k) + +, Ym(k); k k= 1, 2, N; mat 1, 2, 
in which the coefficients b,,,, 4,, K, N, M can estimated. 


= 


um of a periodic 
ment over a 24-h ? 
7 
_ An alternative approach to splitting the demand into a periodic and indepen- § 
_ dently propagated residual component is to consider the total demand L(k) ff 
&§ 
— 


_ methods which can be nied for accurate and computationally efficient forecasts. “4 7 . 
The former is primarily a time-series model to which a wide range of or) 
_ which are well known and well tested can be applied. The Kalman theory = 


is one such approach which can be | extended to simultaneously estimate and 

update time series model parameters, while generating the forecast based on 

the latest state observations (8). The second approach is based on spectral 

_ expansion (12), and decides on the basis of historic data what the most dominant. 

characteristic | modes of demand are. While the we oe for using advanced time 


2 TABLE 1 —Kalman Filter Prediction Model Parameters 


P 


10.9 


4 
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series is limited, this has certain desirable properties 
In the next section, the forecasts of water demands using both of these 
3 techniques will be developed for a given on of water demand data taken from 
the East Worcestershire Water Company (EWWC) in the United Kingdom. oA 
This network has a most advanced computerized telecontrol scheme which has - 
a been extensively described elsewhere (5). The calculations described here have 
i 4 formed the basis for an on-line forecasting scheme which has been developed — 
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"For in an line environment, it is easy to show, (8) that 
isd model is ‘split into independent periodic and residual components, the Kalman 
Bn cg. theory can be used. The demand forecast one period ahead Lk + 

1) is related to an estimate of the state x(k) of a linear dynamical system — 


‘driven by white noise. The state est estimate is such that the error covariance 
matrix 


7 
. 


‘is minimized. The very well known Kalman filter ¢ equations have been documented — 


of 


TEMPERATURE DEVIATION 
FROM NORMAL 


FIG. 2.—Kalman Filter Prediction 
elsewhere (11). The s ante ‘ectimate ¥(k) can be viewed essentially as a measure x 
of the residual component of the demand. It can therefore be used as a means 7 


of updating pumping schedules to account for demand variations about the 


_ Table 1 gives the numerical values for a model of six harmonics with n 
24 =2, m= po! 1, identified using the Flectcher-Powell algorithm (8) for minimizing 
the least squares fit of £(k) to L£(k) for the data shown in Fig. 2. The Prediction — 

accuracy over this period is shown, and the steady state error variance is of — 
the order of 3% of peak consumption. The numerical values ad the cocfiic ficients © 


3 
@ 1 
Time of Prediction: 1:00a.m. Tuesday,August 7th 
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from historic data, it is a relatively easy matter to update this model to take 
account of the latest observations in calculating the latest prediction, €.8., the ; 


an IBM 370/ 165 « computer, while the | generation and updating of forecasts took — 
only 15 sec on the same machine. Another reason for its suitability i is the —_ 


of the noise model ensure asymptotically stable and aperiodic forecasts. Fig. : 
shows som 
_ the model par 
 implementatic 
= a 
| 
= 


WATER SUPPLY NETWORKS. 
that the optimal p pump ‘schedule can be related to the optimal demand estimate — 
&(k) which is fed back from the network. This computation involves the solution 
_ of a matrix Riccatti equation for a linear model-quadratic cost pump schedule 
a which can also be quite modest. Thus, the state estimate ¥ (k) used for forecasting 
4 also be used directly for optimal control. ine 
4 _ By way of comparison, a spectral expansion forecast was carried out | on = 


‘the same data. A record of 53 days excluding weekends and Mondays was — 
used. An exponential smoothing model with constant 0.14 was used to eliminate -~ 
- the base trend in the data, and the une expansion procedure carried out 
_ The accuracy of this method is particularly good for short term predictions 
“up to a few hours ahead. ‘Fig. 4 shows, however, that for 24 h predictions 


there seems to be little to choose between the methods in terms of accuracy. 
_ In terms of its computational requirements, the spectral expansion procedure — 


9 sec CPU time on the IBM which is modest. . Also, 


nw 


« 


> 


4 of P Prediction Accuracies 


> * 


itd does not ; require | any data other than past values of water demand. However, -_ 
a substantial amount of storarage is required for processing these historic data. 
In addition, a completely new calculation must be carried out when new 


__ observations are made, and there is no simple w way of updating | on-line pump > 


controls based on latest observations. For this reason alone, the Kalman filter 
approach i is preferable for on-line 


_ This paper has taken a preliminary | look at the requirements of a a forecasting 
system for on-line pump scheduling in ; a predictive environment on water supply 
networks. Two methods are suggested as being appropriate in terms of ace J 
curacy—the Kalman filter approach and the Karhunen-Loéeve spectral expansion — 
method. The former ® more attractive in terms of on-line implementation a 
the reason. . Part of the forecast involves tl the estimation 


g 
‘Time ahead of prediction ( 
4 


of the state vector of residuals” demand. This gives the deviation of 


actual demand from the periodic base component, and can therefore be used 
on- -line for calculating the optimal pump control. sss 
_ The conclusions drawn in this paper are based on actual data from an advanced = 
_telecontrolled water supply network in the United Kingdom. A forecasting scheme _ “4 


_ based on the methodology outlined in this paper has currently been successfully _ 


implemented on that network, 
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Aprenoix 1_—Notation | 


eee symbols are used in this p 


= Tesidual model 
spectral expansion coefficients; 
model noise coefficients; 
Lk) Water demand in period 
= number of sample in in ensemble; 
= index of sample functions; == 
of days demand; 
number of periodic terms in Fourier s seri 


‘State estimate error covariance; 


= dt covariance matrix; 
Im “mth sample function for water demand; 
deviation from average; 
(k) state estimate of residual 
= periodic component of water demand; 
y, (k) residual component of water demand; demand; 
(k) = spectral expansion error; = 
= ‘ith eigenvalue c of Q; and 
(k) = ith ‘eigenvector of Q. 
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~CONTAMINANT TRANSPORT Move, a 


Yasuo Onishi," M. “ASCE vit 


‘The effects of man’s activities on surface- water systems must be anticipated | 


© prevent undesirable environmental consequences. Sediment movement is an 
particularly important mechanism that must be assessed because sediment is 

a major carrier of many hazardous materials. 

_ Mathematical modeling is a valuable method for predicting sediment response © 
and contaminant migration and fate. However, most available models do not 
account for sediment / contaminant interactions (I, 10,18). ‘The ‘important interac- 
tions include contaminant adsorption by and « desorption from sediments, and 
the migration (transport, deposition, scouring) of contaminated sediments (2,6, 


Bt Mego when these mechanisms are considered can accurate predictions 
a This paper considers the formulation and the verification of the sediment 


and ‘contaminant transport 1 model, _FETRA. Because _FETRA accounts for 
sediment /contaminant interactions 1), it provides n more complete. information 


than previous models. FETRA was applied to the James River estuary, Virginia, i 
to simulate sediment movement and the transport of the pesticide Kepone which : 


was discharged to the river in substantial pened during the early 1970s (17). . 


FETRAi is an two- ‘that utilizes the finite e clement 


computation method with the Galerkin-weighted residual technique. The following 


three submodels were coupled to account for the sediment / contaminant interac- 
tions: A ‘Sediment submodel, (2) a dissolved contaminant transport 


“contaminants a are those adsorbed by sediments. 
_ Sediment Transport Submodel.—Important sediment characteristics such as 


“fall velocity, critical shear stresses for erosion and deposition, and adsorption — 


capacity vary significantly with sediment sizes and types. Accordingly, sediment — 


a "Staff Engr., Pacific Northwest Lab., Battelle Memorial Institute for the Dept. of Energy, 


P.O. Box 999, Richland, Wash. 
aH - Note.—Discussion open until February 1, 1982. To extend the c closing date one month, 
a written request must be filed with the Manager of Technical and Professional Publications, — 
ASCE. Manuscript was submitted for review for possible publication on June 11, 1979. 
This paper is part of the Journal of the Hydraulics Division, Proceedings of the American 
Society of Civil Engineers, © ASCE, Vol. 107, No. HY9, September, 
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movements and contaminant transport are separately for 
each sediment size fraction or sediment type. (Currently, FETRA handles three 

- sediment size fractions or sediment types.) The sediment transport submodel 
_ includes the mechanisms of: (1) Advection and diffusion / dispersion of sediments; 

(2) fall velocity and cohesiveness; (3) deposition on the riverbed; (4) — 
the riverbed (bed erosion and armoring); and (5) sediment contributions 
point/nonpoint sources, and subsequent mixing. Sediment mineralogy 

—_ quality effects are implicitly included in Items 2, 3, and 4. This submodel j 
also calculates the changes in pre riverbed elevation 


“The governing g equation n of the averaged, ‘two- dimensi en 
~ transport for jth sediment size fraction or sediment type was obtained joel 
following three-dimensional mass conservation equation (3, 16): over OF 


oy oz / 


in which G= = concentration of sediment ‘of jth” size fraction, or type (weight 
e of sediment per unit volume of water); Q., = source strength of jth sediment 
- contribution; ¢ = time; U = velocity component of longitudinal (x) direction; 


V = velocity component of lateral (y) direction; _W = velocity component — 
of vertical (z) direction; fall velocity of ‘sediment particle of jth size 

fraction or type; x,y,z = = longitudinal, lateral, and vertical directions i in Cartesian 
coordinates, respectively; and ¢,,¢€,,€, = diffusion coefficients « of longitudinal, 

lateral, and vertical directions. The boundar conditions were _ Ae 


z=0 .. 


+ €, 


ay C= =C, at y= 0 and B. 


in which B = width of the river; yo = constant concentration of jth sediment; 


he = “depth; Sy = sediment deposition | rate per unit bed surface area for jth 
d sediment size or type; S,, = sediment erosion rate per unit bed surface —_— ; 


/ 


for jth sediment size or type; = coefficient, i.e., probability that particle 


settling to the bed is re- sansaiied: ; and Gy = lateral influx of jth sediment. 
For the model development work y was assumed to be unity, i.e., for the 
- same . flow condition, all suspended matter ‘settling on the riverbed as a result 
of fall velocity alone returns to the flow. The vertical flow velocity, W, was 


Eq. | was vertically integrated to obtain the two-dimensional sediment- -transport — 


equations. ‘Since velocity and sediment are generally 
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not uniform entdioiine: we use the following. approach sim similar to to the one e used 


and lateral velocity, respectively; and = fluctua- 
_ tions from the depth averaged values of concentration of jth sediment, longitudinal r 
velocity, and lateral velocity, respectively. Note that c”,u”,v” are spatial — 
deviations, not temporal deviations as are usual in turbulence ; analysis; all of 
the temporal averaging has been carried out prior t to writing Eq. :. —_ 
As in the Boussinesq diffusion coefficient concept let 


= yh = —hD 


in which D, and D, = the dispersion coefficients of x and y y directions. 
_ The dispersion coefficients of x and y directions were assumed to be the 
Z “same for all sediments and contaminants. Noting the equation of continuity, 

the kinetic water surface boundary condition, and Eqs. 2, 3, 9, and 10, the 
“ following final expression of sediment transport was obtained by substituting, 

the ee een into Eq. 1 and by integrating Eq. 1 over the 


+ 
n which = D.; and D,. The finite element 


method was to solve Eqs. 4 and 


one To obtain the sediment concentration, sediment erosion and deposition 
1 _ ‘fates, S,,and Sp, ,» were calculated as described as follows. At each nodal point, Z 
Sy and Sp , were evaluated separately for each sediment size fraction or sediment _ 
type becouse erosion and deposition characteristics are different for cohesive 
_ and noncohesive sediments. Both erosion and deposition of noncohesive sedi- 
7 ments are affected by the amount of sediment the flow is capable of carrying, Z 
_e.g., if the amount of sand being transported is less than the flow can carry 
for given hydrodynamic conditions, the flow will scour sediment from the stream a 
bed to increase the sediment transport rate. This process occurs ‘until the actual 
ecomes equal to the carrying capacity of the flow 


= 
: 
. 
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or or until all the available bed sediments are scoured, whichevea is first. dete 
the flow deposits sand if its actual sediment transport rate is above its capacity 
_ to carry sediment. Because of the simplicity of the formulation, DuBoy’s formula 


(20) was used to the total sediment ‘transport: capacity = 
in which +, = bed shear = critical shear stress defined determined 
by DuBoy (20) as a function of sediment size; and J, = coefficient defined 


and determined by DuBoy (20) as a function of sediment size. Although the i. 
7 DuBoy’ s formula is sometimes classified as a bed-load formula, it has been : 


widely used to. calculate. the total sediment load because of its simplicity (20). 


_ Comparisons of computed and measured sediment loads for the Niobrara * el 
r near Cody, Neb. and the Colorado River at Taylor’s Ferry, Colo., reveal that 
- DuBoy’ s formula tends to overestimate a total sediment load under low flows, — 


but provides better estimates under high flows (20). 
‘The sediment capacity, Or, per | unit at down-~ current: 


may 


Van 


J 


in which AL = diesnen | between the up-current and down-current eae 
and Q,, = is obtained as a product of sand concentration and a flow rate 


TA 


_ For cohesive sediments, e.g., clay and silt, the following formulas developed © 


ba 
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in which | M, = = erodibility coefficient for jth "sediment; = critical shear 
stress for sediment deposition for jth sediment; and Ter = critical shear stress — 
_ for sediment erosion for jth sediment. Values of M,, t.p,, and t,,, must be 
_ determined by field or laboratory tests, or both (8,15, 20). Similar to erosion - 
of sand, cohesive sediment will be eroded with the rate of Sp OF until all 
= the available cohesive sediment in “the bed is scoured, — results in 
When the fall velocity, W, (as shown 
= and 1 no aggregation occurs, the fall could be assumed 
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TRANSPORT MODEL 
in which constant depending on the sediment type (8). 
i Some : studies (2,17) reveal that because of its large adsorption capacity, organic — 

_ matter is a very important carrier of contaminants. Unfortunately, there have 
not been enough studies to quantify the rates of transport, deposition, and 
erosion of organic materials when these materials are transported other than 

_ by attaching to cohesive sediments. Since the mechanisms governing the erosion 

deposition of organic "matter somewhat similar to those for cohesive 
sediments, Eqs. 15 and 16 were also used to solve the erosion and deposition 
rates of organic matter. The selection of the values of W,, M,, t.p,, and T.,, 
in Eqs. 15 and 16 should reflect the characteristics of these materials in terms 
their density, size, cohesiveness, and consolidation. 
- To simulate riverbed conditions as regards the bed’s elevation change, sediment 

- Gintiation, and armoring, FETRA divides the riverbed—except for the top 

_ layer—into a number of layers with a standard thickness. The thickness of — 4 

the top layer is equal to or less than the standard thickness at any given time. 

Each layer consists of a combination of clean or contaminated sediments of 
three sizes (or types), or both, selected for sediment transport as mentioned | 
‘previously. Based on the Sediment erosion or deposition rate Sx or 
calculated | by Eqs. 13 ‘through 16), ‘sediments of each size fraction (or type) 
will be scoured from or deposited to the bed, changing the thickness of the — 

_— top layer and possibly the number of bed layers. When the top bed layer thickness ; 
_ reaches more than the standard thickness due to sediment deposition, a new 

top layer will be formed. On the other hand, when all the sediment in the | 

bed layer is scoured, the bed layer immediately below the: original 
layer becomes the new top layer. This process will continue until the actual — 
sediment transport rate becomes equal to the carrying capacity of the flow 
or until the available bed sediment is completely scoured » whichever occurs — 
- For: Gen sediments of the second and third size fractions (or types), ‘the number 7 
of bed layers eroded cannot exceed the number of layers eroded for the sediment — 
os of the first size fraction (or type). In other words, sediments of the first sediment 
size fraction (or type) cover and protect the sediments of the second and the © 


within the were also” obtained by keeping track of ‘the amount of 
period due to erosion or r deposition of contaminated sediment, or both, 
_ direct adsorption / desorption between dissolved contaminants and bed sediments. 
_ Dissolved Contaminant Transport Submodel.—Dissolved contaminants interact © 
- both with sediments in motion (suspended and bed-load sediments) and with 
Stationary sediments in the riverbed. To account for these interactions, this" 
submodel includes the mechanisms of: ray) Advection and diffusion /dispersion 
of dissolved contaminants; (2) adsorption (uptake) of dissolved contaminants 
+: both moving and stationary sediments or desorption from the "sediments 
into water; (3) chemical and biological degradation, or radionuclide decay of 
and (4) contaminant contribution from point/nonpoint sources, 
subsequent from wastewater _Giacharges,  overtend 
runoff flow, fallout, 
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a part € point/nonpoint source Contributions. ects Of wa er quality, 


e.g., pH, water temperature, and ‘salinity, and sediment characteristics, e. g. 7 
clay minerals, can be included by changing the distribution (or partition) 
governing equation of dissolved contaminant 1 transport for ‘the: three- 


ax 


€ K,(C, K4G, 
In addition to the previously defined symbols: K,, = distribution (or partition) 
coefficient between dissolved contaminant and particulate contaminant associated 
with jth sediment; K,= rate of f contaminant: adsorption or desorption to reach 
an ‘equilibrium « condition with jth sediment in motion; G,= - particulate contaminant _ 
concentration associated with/th sediment (weight of contaminant or radioactivity 
per unit volume of water); G,, = dissolved contaminant concentration (weight 7 
of contaminant or radioactivity per unit volume of water); Q,, = source strength 
of dissolved contaminant; and A = chemical and biological degradation rates, 


or radionuclide decay rate of contaminant. The distribution coefficient, 


K 


in = ‘fraction of contaminant sorbed by jth fraction 

; of contaminant left in solution; M; = weight of jth sediment; V,, = volume 
of water; and G,/G,,. Eq. 19 may be rewritten as 
of Adsorption of contaminants by ‘sediments or desorption from ‘sediments a are 
assumed to occur toward an equilibrium condition with the rate, K, (unit of | 
recriprocal of time), if the particulate and dissolved contaminant conesntentions 
differ from their equilibrium values as expressed in Eq. 20. The rept 
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at 


in which D, einen diameter of particulate 
‘contaminant ‘concentration associated with the jth sediment in the riverbed (weight _— 


| 
act 
q 
d 
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q = 


i Lye 
of contaminant or sais per unit weight of sediment); K,, = rate of 
Es contaminant adsorption or desorption to reach an equilibrium condition ' with = 


ee o~ riverbed sediment; GC wo = constant concentration of dissolved contaminant; 


~" seaialaane and stationary jth sediment in the riverbed, or both. The distribution 
coefficients for both moving and stationary sediments were assumed to be the 
; ay as long as these sediments have the same characteristics, e.g. , diameters, — 
_ Clay minerals , organic « content. The difference in the adsorption /desorption 
ee - mechanism between these two yo sediments is the time required to reach equilibrium. = 


We also assumed that direct adsorption or desorption, or both, from bed aa aa 


occurs only with sediments on the bed surface. 
order to integrate Eq. 18 | 


= 


_ in which G, = depth averaged value of particulate contaminant concentennen 

associated with jth sediment; G7 = fluctuation from the depth averaged value 
of particulate contaminant concentration; 6. = depth averaged value of ¢ dissolved a 

contaminant concentration; and = from the depth averaged value 7 


equation of dissolved contaminants: 
0 


- K, 


ay ax 


we 


ay. 


K,K4C,+ > Kay, 


boundary for this e are the same as those in Eq. 23. 
Particulate Contaminant Transport Submodel.—The transport of 

- attached contaminants is solved separately according to sediment size fraction — 
types). submodel also includes the mechanisms of: (1) Advection and 
diffusion / dispersion of particulate contaminants; (2) adsorption (uptake) of 
dissolved contaminants by sediment or desorption from sediment into water; = 
(3) chemical and biological degradation, or radionuclide decay of contaminants; ‘ 
, +s (4) deposition of particulate contaminants on the riverbed or erosion from the ‘y 
riverbed; and (5) contaminant contribution from point/nonpoint sources and © 
ixing. Again, contributions to a river system f from wastewater 
discharges, overland runoff flow, fallout, and ground water» may be treated 
_ as a part of the point/nonpoint sources contributions. The three- dimensional — 
distribution of the particulate contaminant within the riverbed is also 


and = lateral influx of dissolved contaminant. 
| 
E- 
= +—(K P)— 
wee 
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"sional transport equation for contaminants by the jth sediment may 


8G, 4G, 
in which Q, = source ce strength c of pli ae contaminant associated with jth 
sediment. The boundary conditions for this case are 


_in which q, = lateral influx of particulate contaminant associated with jth sedi- 
ment; and Gi 


of the particulate contaminant transport equation: 


The a conditions for "this case are the: same as those expressed in Eq. 

’ * Finite Element Method. —Because se of its increased solution accuracy and ready _ 
accommodation to complex boundary geometry, the finite element solution BP: 

: technique with the Galerkin weighted residual method (4) was used in this study 
to solve Eqs. 11, 26, and 31 with the boundary c conditions of Eqs. 4, , 23, and 

_ 30. The flow domain was divided into a series of triangularelements interconnected = 

: at nodal points. Six nodes were associated with each triangle, three at the 


vertices and three on the midsides. A quadratic approximation was made for 7“ 
the sediment and contaminant concentrations with each element. Either linear 
or quadratic interpolations can be used for variations of flow depth and velocity 
_ within an element. The computer program was written in FORTRAN IV language > 
to adapt the model to CDC-7600 and VAX computers. During the following 
~ model application, the computer time required to calculate all seven substances 
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computational node per time : step was 0. 0028 cp the CDC-7600 


Move Appucation > | or be 


The FETRA code was ‘to the of sediment and 


Kepone, a pesticide released to the James River estuary (Virginia) through Bailey / 
- Crock as shown in Fig. 1 (2,7,17). An 86-km (53- ~mile) estuarine reach between © 
City Point (River Kilometer 123 or River Mile 76) and Burwell Bay (Rive 


37 or River Mile 23) was simulated. Tidally influenced depth and 


velocity distributions in the study a area were obtained by the unsteady, one-dimen- | 
p _ sional code, EXPLORE (1), which gives cross-sectionally averaged velocities a 
and flow depths. These one-dimensional results were utilized by FETRA to 
obtain distributions of sediment and Kepone. Thus, results presented 


BAILEY CR q 


be si joomibet b > bak 


bun sin) FIG. 1.—James River Estuary 
_ Applying FETRA to the estuary also involved model calibration and verification b 
Since field data (2,17) indicated that a considerable amount of particulate Kepone 
_ was transported by organic materials moving independently with other sediments, 
sediment transport was modeled for three sediment types: (1) Cohesive sediments 
(clay and silt); (2) noncohesive sediments (sand); and (3) organic matter. The -, 


modeling procedure at each point in time involved first simulating the transport, q 


‘deposition, and erosion of sediment. The sediment transport results were then 


J ‘used to simulate dissolved and particulate contaminants by accounting for the - ’ 
= contaminant interactions with the sediments. Changes in riverbed conditions 
_ were also calculated in terms of: (1) Changes in the river bottom’s elevation; 
ie (2) the distribution of cohesive sediment, noncohesive sediment, and organic 

- matter in the bed; and (3) the distribution of pasticuinte Kepone in in a the bed. wer 


a 


| 


= were abunet to fit the simulation results to the measured data were the © 


dispersion coefficients, plus the erodibility coefficient and critical shear stresses 


used to calculate deposition and erosion rates of cohesive sediment and organic a 
matter (Eqs. 15 and 16). The other parameters, e.g., Kepone distribution 
Coefficients sediment sizes, sediment fall velocities, and the like, were determined © 
by theoretical and experimental analyses, and field measurements prior to the 
| simulation. Thus, the “major” calibration effort _Teproduced sediment 
"distribution patterns similar to the measured longitudinal distribution of sediment 
concentrations for the 86-km (53-mile) river reach (2,11) by adjusting © the 
longitudinal dispersion coefficient, erodibility coefficient, and critical shear 
stresses. The lateral dispersion coefficient was assumed to be zero because 


Ideally, FETRA should be calibrated with two- dimensional data. However, 
the use of the one- -dimensional, hydrodynamic ‘model, EXPLORE, "precludes ~ 
two-dimensional testing. The results mentioned in this } paper are those obtained 
onthe 3lstday ofthe simulation, 

| _ The model parameters and coefficients were adjusted until the model repro- — 


duced measured total sediment concentration (the sum of cohesive sediment, 
organic matter and ‘sand | being transported as suspended and bed loads) at 
maximum ebb, slack tide, and maximum flood for a net freshwater river discharge _ 

of 58.3 m*/s (2,050 cfs) at City Point (River Kilometer 123 or River Mile 76). rf 

& 2 shows changes of computed sediment concentrations with time at River 7 
Kilometer 75.7 (River Mile 47. 0). In this figure, the maximum ebb occurs at | 

__ approximately 3 h and 15 h. The maximum flood occurs at approximately 8 : 

hand 21h. Fig. 3 shows tidally averaged, longitudinal distributions of computed 
sediment concentrations for each sediment component (cohesive sediment, 
organic matter, and sand) and total sediment, together with measured total 
sediment concentrations (2,7). The longitudinal variation of sediment concentra- 

tion is relatively small due to the lack of a clear null zone (where fresh and © 
Saline water intermix). As shown in Fig. = ‘the agreement between the computed 

~. erification.—The sediment transport submodel was verified for a net freshwa- 4 a 


_ 


ter discharge of 247 m’/s (8,700 cfs) without changing the model parameters 
_ and coefficients. Computed sediment concentrations were compared with mea- 
sured total sediment concentrations for various tidal stages. Slack tide and 
tidally-averaged computer results $ are shown in 1 Figs. 4 and 5 together with field BS 
data (9). These figures include computed sediment concentrations of each type 4 
of sediment (cohesive sediment, organic matter, and sand) and total sediment 
(the sum of the sediment components). 
_ Computer results revealed that there are two peaks of the total sediment 
concentrations of approximately 100 mg / One occurs between River Kilometers 
75 and 95 (River Miles 47 and 59), and the other around River Kilometer 55 


(River Mile 34). The first peak was produced by extensive river bed scouring o 


a 


in the vicinity of this reach, and the second was due to the existence of a 
clear null zone. In Figs. 4 and 5 the computed results are compared with field 
data obtained in 1965 (9) and 1966. The 1966 data were supplied by M. M. 
Nichols, Virginia Institute of Marine Sciences (vane. As round by Figs. 4 
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To test the contaminant submodels, ‘computed results for Kepone © 


compared with measured data at various tidal stages (2,7). Because the 
= a s 

TOTAL SEDIMENT | 


~~ — COHESIVE SEDIMENT wer aft 


+ 


/ 
2 2 
: 
FIG. 2. —Changes of Computed Sediment Concentrations with Time at River KRometer 


75. 7 (River Mile 47.0) for Net Freshwater Discharge of 58.3 m 3/s (2,050 cfs) _ wie 


- ORGANIC MATTER 


= 


s & 


(DAWSON, 1978) 


an FIELD DATA (HUGGETT, 1978) 
4 


ENT CONCENTRATI 


a DALLY AVERAGED SEDIM 
2 


A 

| 

KILOMETERS 


FIG. 3. Tideity y Averaged Sediment Consenwation of Each Sediment 
for Net Freshwater Discharge of 58.3 m*/s (2,050 cfs) wel 
i adjustable model parameter and coefficient are all related to the sediment transport 
a] as mentioned are no parameters and coefficients 
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‘Thus, no was performed for these contaminant transport si submodels. 

Predicted particulate Kepone concentrations associated with each type of 

| Fred and weighted average particulate Kepone at the maximum ebb tide a 

_ Shown together with measured data of average particulate Kepone concentrations 
Fig. 6). The computed results and the field data closely agree. its 

5 _ Most of the dissolved concentrations that occurred in the reach of the James 

River under study were below the the detection limit of approximately 0. .005 pe/L 
= 0.010 ug/L (2,7). From the measured particulate Kepone concentration and ae a ; 


the value of the distribution coefficient, the expected dissolved Kepone con- 
centration is somewhat below the detection limit but is roughly the same patel 


4, 


ngitudinal Distributions of Computed Total, Dissolved, and 
Concentrations Maximum Flood “ai Net Freshwater of 4 


~ computed dissolved, particulate, and total (sum of dissolved and particulate) ony 
-Kepone concentrations with time at River Kilometer 75.7 (River Mile 47.0) | 
are shown in Fig. . | Fig. 8 shows the computed dissolved Kepone concentration — 
_ together with particulate and total Kepone concentrations at the maximum flood 
tide. This figure indicates that dissolved Kepone concentrations vary from 
approximately 0.0048 yg/L-0.0084 we/L. concentration of dissolved 
Kepone predicted by the FETRA code is the highest possible value found in an : 
the study area, but is still below the detection limit. Total Kepone concentrations a 
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atin of Predicted Riverbed Elevation Changes Due to Sediment Deposi- 
tion or Bed Scouring, or Both, for Net Freshwater Discharge of 58.3 m */s (2, adi 
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“Kilometer 77 Mile 48) at t the flood tide (Fig. 8). ‘This figure 
also suggests that the suspended sediments carry 12%-53% of the total Kepone 
being transported, while 47%-88% is in a dissolved form. The mathematical 
atten results also show that at Burwell Bay (River Kilometer 37 or River 
Mile 23), the total Kepone concentration averaged over a tidal cycle was, $s an 
estimated 0.0076 g/L. This concentration implies that 14.0 kg of the Kepone 
_is transported annually from the tidal James River estuary toward Cheseapeake © 
4 Bay for the net freshwater discharge of 58.3 m’/s (2,050 cfs). Of this total, 
2.2 kg of the Kepone is carried by ‘Suspended sediment, while 11.8 kg of the 
Kepone migrates from Burwell Bay e each year in a dissolved form nee this 
_ Two computation examples for riverbed conditions are shown in Figs. 9 and ie 
10. Fig. 9 indicates accumulated bed elevation changes predicted after the 
: one- -month simulation. In this figure, positive values along the vertical axis 
indicate the amount of net sediment deposition in millimeters during the one- -month — 
| sition, while negative values show the amount of riverbed scoured. The 


_ geometry of the James River estuary, in 1 which many bays a are € connected by 4 
_ ‘arrow channels. For this flow case, an annual net bed deposition rate averaged 
(over the 86-km (53-mile) reach was predicted to be 1.7 mm. Fig. 10 shows 

” the predicted change in bed-surface Kepone concentrations that occurred during © 

- the one- -month simulation. Definitely, the Kepone level tended 1 to o decrease n near i 


The unsteady two- dimensional ‘finite element model, FETRA, was developed Bt 4 


4 


by including the of sediment /contaminant interaction, | e. g., absorp 
tion and desorption of contaminants; and transport, deposition, and resuspension 
of contaminated sediments. The FETRA code was applied to the James River — 
estuary to simulate the transport of sediment and the pesticide, Kepone. Although , 

_ additional field data were needed to rigorously verify the model, predicted © 
concentrations agreed well with field data obtained i in the ee 

ies 
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(DOE), whose Pacific Northwest Laboratory (PNL) is operated by = 
Memorial Institute under Contract DE-AC06-76RLO 1830. . The. application of 
the model to the James River estuary was sponsored by the» United States - 

Environmental Protection Agency (EPA), through an interagency agreement 
_ between EPA and DOE. The writer wishes to express his gratitude to R. -. 
_ Baca for his assisting in the model development : and to S. E. Wise for sameins 
in the model 
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The following symbols ai are used in this nee: 

= concentration of sediment of pe t pe (weight of sediment 


unit volume of water); 


depth averaged concentration of sediment for jth type; 


= fluctuations from depth concentration of of 
constant concentration of jth sediment 
= dispersion coefficients of x and y directions; 
, = fraction of contaminant sorbed by jth sediment; — Snees* 
fraction of contaminant left in solution; 


in river bed (w weight of contaminant or radioactivity p per unit weight 


particulate contaminant concentration associated with jth sediment 


depth averaged value of particulate contaminant concentration 


radioactivity per unit volume of water); ren 
depth averaged value of dissolved contaminant concentration; 
” = fluctuation from depth averaged value of dissolved ‘contaminant 


ak 


= constant concentration of dissolved contaminant; etree 
= flow depth; 
of contaminant adsorption or desorption to reach an 
um condition with jth riverbed sediment; 
= distribution | (or partition) coefficient between dissolved contami- = 
| nant and particulate contaminant associated with jth sediment; a 
Aa K, = ‘fate of contaminant adsorption or desorption to reach equilibrium _ 


coefficient for se sediment | of jth size fraction or type; 


actual total sediment load per unit width; 
strength of particulate contribution ass associated 


Der 

= 

| 

— 
| 


Q, = total codiment transport « capacity of flow per unit width; 


source strength of dissolved contaminant contribution; 


lateral influx of the jth sediment; so 7 
influx of dissolved contaminant; 


lateral influx 0 of contaminant associated with | jth sedi- 
siento deposition rate per unit bed surface area for jth sediment; 
sediment erosion rate per unit bed s surface : area vines sediment; 


component lateral (y) direction; 


depth averaged lateral velocity; 
lateral velocity fluctuation from depth lateral v 
_ velocity component of vertical (z) direction; 


P s longitudinal, lateral, and vertical directions in Cartesian coordia- 
= coefficient, i.e., aed that particle settling to bed is re-e -en- 
density of jth sediment in river bed; 
distance between up-current and down- current locations; 


= diffusion coefficients of lateral, and vertical diree- 
rate of contaminant; 


= 
bed shear stress; 


critical shear stress for sand defined by DuBoy; a 
= critical Shear stress. for deposition for jth cohesive 


itical shear stress for sediment erosion for jth cohesive sedimen’ 


TRANSPORT MODEL __ 1107s 
= 
longitudinal velocity fluctuation from depth averaged longitudinal 
Ww 
> 
‘ 


= 
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a 
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idth ratio of 1.0, whick were leben 
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: - applies separately to tables and references. In referring to a figure, table, or reference that i 
appeared in the original paper /technical note use the same number used in the original. = 
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Publications of A ASCE fo for ‘more detailed information on preparation and submission of 
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(CALCULATION OF STRONGLY CurveD OPEN F Low 


Closure by Michael Leschziner® and Wolfgang Rodi PRED! 


‘The v writers ‘like ‘to ery Gotz, for t their i interest 
‘in the paper ai and for their comments. Gotz has pointed « out that curved open 
channel flow with a ratio of radius-to-width <3.5 cannot exist without large 
separation zones. However, the measurements of Rozovskii (18) for a radius-to- _ 
_ idth ratio of 1.0, which were taken for validation of the mathematical model — 
in the paper, do not show any reverse-flow regions, i.e., where the longitudinal 
"velocity is in the upstream direction. The occurrence of flow separation depends 
also on parameters other than ‘the: radius- -to-width ratio, and the “experimental 
flow situation of Rozovskii was chosen for the model validation precisely because 


The main comment of Anwar and Rodger is that the empirical coefficient, 
% in the turbulence model is not really a constant as assumed in the paper. _ 
The writers are well aware of this fact and should like to point out Refs. 42 
and 43 where refined turbulence models are described that simulate the ae | 
s c,, near the free surface which is necessary in order to obtain the well-known 
"parabolic eddy-viscosity distribution in open channel flow. However, even though © 
c,, can vary significantly across the flow, and its average value also somewhat _ 
; from one flow to another, in many cases the calculations are not very sensitive _ 
to the c, distribution so that the main features of a great many flows can in 
fact be predicted well with a constant value of é = 0.09. This is demonstrated _ 
in Ref. 44 where results obtained with the k- -€ model are compared with 
experimental data for many different hydraulic flow situations, thus refuting — 
the discussers’ statement that the k-e model cannot be applied to flow probl2 >ms- 
more general than isotropic turbulence without tuning the constants. The fixed 
constants of the k-e model are, of course, not suitable for each and every © 
_ ‘ flow, and exceptions ; such as strongly buoyant flows are also described in Ref. 
: 44. This reference also concludes that simpler models like the Prandtl mixing 
4 length theory work fairly well for simple channel flow and submerged jets bet 
require different empirical input for different problems and are not suitable — 
7 whenever turbulence transport or history effects are important and when the 
“flow i is more complex than a shear layer, e.g., ‘recirculating f flows. In such 
2 cases, the additional computational effort required to solve the k and € equations, 
which is less than 20% with a suitable numerical scheme, is certainly justified. 
Calculations of strongly curved open channel flow as described in the paper 
7 are not particularly sensitive to the turbulence model employed because in andl 


case the distribution of the velocity components, on which attention was focused, 
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“October, 1979, by Michael Leschziner and Wolfgang Rodi (Proc. Paper 14927). renal 
P.O. Box No. 88, Manchester, England, M601QD. 


ag Research Engr., Univ. of Karlsruhe, Karlsruhe, W. ‘Germany. 


ry 
| 
| 
7 


governed mainly by pressure forces. Thus, a simpler model would 
probably have yielded similar results, including the vorticity model of Prandtl 

(41). The k-e model was | used simply because it is at present one © of the most — 

Anwar and Rodger have based their comment on the variation of c , on the 

j evaluation of c, from measurements reported in Refs. 35 and 36. According - 

to Eq. 6, evaluation of c , Tequires knowledge of the eddy viscosity v,, the | 
turbulent kinetic energy k, and the dissipation rate «. There is no indication 
a. Ref. 35 that « has been measured, and it is therefore not clear to ) the writers if 
& the discussers have actually evaluated c,. This casts some doubt on the 7 


numerical values for c,, quoted by the discussers, but the comment that c,, 


42. Hossain, M. S., ‘‘Mathematical Modelling of Turbulent Buoyant Tie (in German), j 
_ thesis presented to the University of Karlsruhe, at Karlsruhe, Germany, in 1980, 
partial of the degree of Doctor of Engineering. 
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~ _Channel Flow with an Algebraic Stress Turbulence Model,’’ Report SFB 80/T/ 187, &§ 
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44. Rodi, W., Turbulence Models and their Application in Hydraulics—a State of a - 
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APPROXIMATE FLoop Routinc Metuops: A Review" 
by ?. i ) 


The thank Ponce and Tan, Price, and Kundzewicz for their 
‘discussion, which clarifies some of the problems dealt with in the paper and a 
_ The terminology for the methods that neglect the inertia terms of the momentum | 
- equation warrants some further comments. In the paper, the term approximate — 
dynamic model was reserved for methods based on the numerical solution of © 
the system of Eqs. | and 2, but which, for reasons of computational efficiency, 
- As pointed out by Ponce and Theurer, the term kinematic model corrected — 
for dynamic effects, used for the Koussis model, is not strictly correct. This 
_ term would be appropriate for a model based on the full equation of the looped 


rating curve, rather than on the Jones formula incorporated in the Koussis — 
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Expressions for a fully dynamic curve written in the 


_ changing flow rate at one section have been proposed by Fread (28) and Simons, p 
Li, and Simons (29). Because of its inability to model the effects of downstream 
_ disturbances and some of the approximations made in its derivation, a model 
based on such expressions would still be kinematic in nature, though corrected 
for dynamic effects, thus justifying this classification. In this context, a more 7 
appropriate terminology for the Koussis model would be kinematic model corrected — 
for diffusive effects. ip oni! 
The question raised by Ponce and Theurer regarding the generality of the : 
+ finite difference scheme proposed by Koussis has been clarified in a Tecent — 
paper by Koussis (30). With regard to the two alternative expressions for the 
computation of the weighting coefficient 6, the writers’ experience ocnties 
_ the conclusion reached by Koussis that the more complicated expression is . 
- ‘The writers appreciate the comments made | by Ponce and Theurer and by 
_— Kundzewiez on upper and lower limits for reach length Ax, and welcome the 
7 opportunity to add some clarifying remarks to this section of the original paper. — 
As pointed out by Ponce and Theurer, a lower limit on Ax is equivalent — 
a lower limit on 8, ie., => 0. The writers agree with Kundzewicz that 


and, like Ponce and Theurer, the first author has successfully used negative 
values of @ (15) . However, as @ and Ax are related through Eg. 15, or Eq. © 
& it is generally possible to increase Ax, while ‘decreasing the number of — 
subreaches, until a positive value of 6 is obtained, without significantly affecting — 
__ The question of an upper limit | on Ax, or an equivalent lower limit on At 
for given Ax, is of more direct importance for practical routing computations. 
In contrast to the theoretically based expressions recommended by Ponce and 
- Theurer and by Kundzewicz, the writers’ condition for Ax, , Eq. 19, is not intended | 
to suppress negative outflows completely, but merely as a guide to keep these 
within acceptable limits. The reach lengths computed from Eq. 19, therefore, 
Bs should be considered as absolute maximum values. The equation also does 
not reflect the interdependence between Ax and At, which is accounted for 
a 7 ww a separate : condition. As stated i in the paper, values of At considerably smaller 
than Ax/(c) =  K may lead to negative outflows. _ The routing results presented — 


in Fig. 6 of Kundzewicz’ discussion appear to have been computed with a 


time step much smaller than his K = 30 min, and the negative outflows ~~ 


| ae thus not unexpected. Also, as Kundzewicz points out, the initial slope of the 
_ inflow hydrograph has an effect on the occurrence of negative ordinates. Our 


~ 


formula was based on experience with hydrographs having the usual S-shaped J 


rising limb. Its slope is reflected in tthe time of riseT,, 42 
comments by Price on a more general form of the diffusion equation 
are a valuable addition to the paper. Similarly to the methods based on the 

; kinematic wave equation, the classification of the various forms of models based 7 

7 “on the diffusion equation depends on the method used for the computation 

Of the routing coefficients c and D. In the form defined by Eqs. 26-29, the | 
= diffusion model is indeed “equivalent to the Koussis model and may, in fact, 


J 
on the weighting coefficient depends on the point of view. For a black box 7 
| 
| 


the | programs appended to 14 and 20. 

8 The writers agree with Price that the concept of an equivalent channel, as 

- introduced by Price and Kawecki (21), is a very useful one for many practical S 

_ applications. There appears to be no reason why this concept could not be © 

_ applied with the Koussis model as well as with other approximate routing methods. 

. The discussion and the references quoted by Kundzewicz draw attention to 

- i another important line of comparative investigation into different flood routing 7™° 

~ methods, based principally on linear systems theory. This approach and the 

_ approach based on numerical modelling techniques, as used by the writers, 
should be seen as complementary and advances with both methods should i help ‘ 
to further rationalize flood routing computations for for practical | applications. 

Errata.—The following corrections , should be made to the original paper: — ; 
Page 1527, paragraph 3, line 8: Should read = instead of “(/e) 
— 1530, line 2: Should seat “ay” instead of **(6) 
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30. Koussis, A. D., ‘“Comparison of Muskingum Method Difference Schemes,” Journal 
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Sreapy-STATE EsTIMATION 0 OF COoLinG Ponp PE RPORMANCE” 


Closure by Gerhard H. ASCE and Masataka Watanabe© 


= 


Adams’ comment t regarding the the a applicability of s of simple steady- State estimation 
= to transient cooling pond conditions is welltaken. 
In general, cooling ponds experience transient effects because of two external 7 
eee The first arises from variable meteorological conditions for which 
the time scale is H(pc, K)~' as mentioned by Adams. The second arises from 
variable power plant operation, and its time scale is the pond residence time - 
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DISCUSSION 


yy QO. in which pond volume, as by the authors. Hence, if 

actual transient cooling pond performance is to be modeled by the sequential 
application of steady-state models—with appropriately averaged or filtered input _ 

; data a as suggested by Adams and Koussis (6)—then, presumably, both time 
scales ought to be considered in the averaging or filtering procedure. In ai 
case, the question remains whether a fully transient simulation (4) is not mere : 


‘Straightforward, more accurate and more efficient. 


‘tem em Bar Bar Resistance oF GraveL-Bep STREAM 


Th chace of Discussion by Arun Kumar* and Rema Devi* 
Fi methodology i is proposed estimating g resistance characteristics ¢ of inactive 
7 gi gravel-bed streams. It is well recognized that, at low depths of flow, lower 7 
ay velocities are observed than the logarithmic law predicts and friction factor — 
_ imereases by as much as 50% (35). However, difficulty arises in estimating 
it quantitatively. ‘The authors’ attempt in this direction is noteworthy. 
‘The authors’ analysis is based on two hypotheses. Observed data that do 
: not fit into the framework of the hypotheses, are rejected, so much so that, 
out of 46 reaches initially chosen, only 15 have been used for further analysis. 
_ Even out of these selected reaches, some points are deleted in order to improve 
correlation between +% and C, (Fig. 6). Rejecting two thirds of the data simply 
_ because of not fitting into the logic of the hypothesis casts doubt on the validity 
of the hypothesis itself. It is contended that field measurement of slope and ; 
grain size distribution are subject to large errors. As these parameters form 
_ part of the model, the resistance can be measured only qualitatively, 
-_ _ Further, the authors do not distinguish between the pool and riffle sections. — 
_ The proposed analysis i is based on the assumption that the flow is approximately : 
uniform. ‘The assumption is valid only att the sections where the 


As 


of the sections. However, at the ‘pool sections, inaccuracies are likely neces 
the flow characteristics are dependent on downstream control. Therefore, the 
writers feel that the large amount of data not conforming to the hypothesis — 
may be partly due to its collection within the pool sections. which We 

The authors have used Keulgan’s Eq. 2 to determine the grain resistance 
with roughness height as 2D,,, and the bar resistance is taken dependent on 
“the grain resistance. Hey (10) has shown that total resistance can be comgund 
directly by taking roughness height as 3.5D,,. Large roughness height is probably 
_ due to development of wake interference losses caused by large — in the ; 


: 
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discharges tally quite pron (standard error in discharges is only +4.7%) 
both for high and low depths of flow. For the illustrative example given by 


- the authors, the ratio of bar resistance C, to the total resistance C at various 
- heights H can be computed by estimating C and C,, from the Keulgan’s equation — 
corresponding to roughness height of 3. SD (D not being known) and 


Except when H = = 0.6 the ratio favorably with authors 


analysis shown in Fig.9(C). 


- Fig. 6 represents a linear functional relationship between log C, and log 


* | . The standard approach of estimating the parameters of linear regression 


assumes that the independent variable 7% can be determined with negligible 
error. However, due to uncertainty in the measurement of slope and grain size 
distribution 1% is subject to large errors. Similarly C, is also subject to large 
errors. In such cases, the regression slope and its confidence limits are to be 
estimated using the approach given by Davis, et al. (36). The writers are of — 
g opinion t that instead of giving a point estimate of discharges at various 
depths, it / would be better, for | given depth slope and grain size distribution, 
if confidence limits of discharges were given. The point estimate gives a false. 
impression that the is highly deterministic, whereas, » real 


36. Bayazit, M., ‘‘Free Surface Flow in a Channel of Large Relative Roughness,” Journal 
Hydraulic Research, Vol. 14, 1976, pp. 115-125. ad 
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i) The authors’ main definition i is based on a stratigraphic distinction (Distinction 4 
“», with reference to the vertical distribution of grain size. The stratigraphic 
definition, although convenient for field application, reflects ideal cases drawn 
from a bed-material composition spectrum controlled partly by provenance. Hence 
the definition may be unduly restrictive or misleading, or both, unless a consistent — 
sampling design is advocated (7, 18), and is not instructive with respect to the 
frequency and manner in | which the protective-layer operates. _exaqel-boronise 
The choice of such process-related, time-dependent terms as ‘‘process of 
armoring,’’ “‘self-stabilizing,”’ and ‘ ‘formation of a stable surface layer,”’ " in 
a many papers, suggests that a process time-based explanation of bed-material 
segregation would be more appropriate for engineering practice. The authors 
continue along such a line, arguing that (Distinction 2) armor and pavement 
could be defined on a basis related to the frequency of disruption of the layer. 
- Unfortunately the authors do not expand on this theme. Presently, such a useful iy 
definition is qualitative and arbitrary, as there are few field data related to “l 
the frequency of disruption of protective surface layers. The authors, neverthe- — 
* less, associate pavement with degraded streams and inactive beds, (22) and © 
a armor with active beds in streams that have achieved some sort of dynamic — 
~ equilibrium (22), a useful concept closely tied to Distinction 2. Isolation of 
: a design flood of specified return period (7) related to bed disruption events | 
‘would be valuable in river stability problems, e.g., bed disrupted by flows > 
-T yr flood = pavement. Beds disrupted by flows < T yr flood = armor. = 


_ [tis not difficult to find an example where the stratigraphic distinction (especially 
‘if "grain “Size composition is } strictly defined) is is ; not in accord with a ‘frequency 


summary remark. An example of a stable cobble bed- surface wherein granules 
Q mm to 4 mm) play an important role in bed st: bility exists (19). The description — 
_ of the bed composition and implied frequency of disruption (19) is not in accord _ 
_—-Bagnold (15) described ‘‘a vaguely defined phenomenon called armoring”’ 


as the continuous erosion of a granular river-bed until the residual surface 
_ solids are so large that further erosion ceases’’—a process which he associated 
specifically with stream reaches below dams. Such a process involves the basic 
elements: (1) Degradation of bed height; (2) spatial concentration of coarse 
particles; and (3) preferential selection by a lag-action of the coarsest clasts 
(17). Bagnold noted that the effect did not apply particularly well to the A meri 
condition of a regular bedload supply through a natural channel. These observa 
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tions suggest that that a number of ' Processes rf result in coarse gravel segregation 

_ The writer (16) suggests that the terms armoring and paving should be restricted 7 
4 _ to active processes that, at some time, have involved the hydraulic entrainment | 
¥ _ and repositioning of a fraction of the coarse bed material such that the sub-surface 

material is protected from entrainment. The individual terms would then be— 
Kote with reference to the frequency of disruption of the layers. Such : a 
definition, would ‘preclude the hydraulic winnowing of fine surface matrix. 
4 materials from around large immobile framework clasts—a process of selective 


- entrainment of certain grain-sizes, the resultant coarse layer being a passive — 
_ response to the process. It is suggested (16) that these latter layers be a 


_ censored-layers. The terminology, being adopted from Tanner (23), encompasses 
the essence of the process involved, i.e., selective winnowing. Using this definition 
the gravel-bed surfaces described by Gessler (3), Harrison (5), and Beschta 4 
& Jackson (14) are not actively formed armor-layers but passively formed — 
- censored-layers. Given certain conditions the matrix material may be replaced © * 
in a censored layer (16,21), strengthening the structure (19), so that no granulo- 
metric distinction can be made between surface and sub-surface deposits. The 
-_ component at the surface will not have altered its function in the over-all 
bed stability, although (using the stratigraphic definitions of the authors) pre- 
viously the bed would have been described as armored and laterally unarmored. 
= importance of being able to identify protective surface-layers, however 
defined, is related, as stated by the authors, to the solution of practical problems 
of river bed stability within a river management framework. The stability problem 
revolves around the accurate prediction of material entrainment and transport. _ 
If a protective-layer controls the release of finer sub-surface materials into 
7 a transporting flow on a rising hydrograph, the progressive disruption of the — 
7 bed should be detected as a flexure or discontinuity in plots of bedload grain-size 
a an entrainment parameter (18)—the discontinuity - defining a critical 
threshold, T yr flood, for protective-layer disruption. It might be expected that 
a threshold for a compact bed surface would be greater than a threshold for 
individual particles (19,20). Nevertheless, the question of bed stability related 
to packing of almost equal size grains remains as crucial and unresolved a 
problem (20) as the stability of bets which develop distinct vertical grain-siz size 
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_ and addressed themselves to the problem of defining what critical experiments 
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at authors have presented an an excellent paper on the use of the terms ‘ ‘armored | 
gravel bed”’ and “‘paved gravel bed.’’ The writer agrees with the usage as proposed 
by the authors because: (1) The suggested usage is reasonably consistent with 
- the usage in other areas of interest to hydraulic engineers and geomorphologists; — 
“and (2) the usage provides information about the physical processes occurring — 
‘There is ene) the use of “armored” and ‘ ‘paved”’ ; e.g., Parker 
= Peterson (14) have used ‘‘paved’’ where the ‘authors used ‘‘armored.”’ In 
discussions the writer has had with others, it became obvious that the = 


. of either “‘armored”’ ” or ‘‘paved’’ could be confusing if one was not quite specific 
as to just what was meant by each term. For this reason an alternative was 


; This suggestion has, in the writer’s opinion, a serious flaw in that it does not 
provide any insight into the physical processes involved, which is one of the 
major benefits obtained from the authors proposed usage. The writer suspects, 
one would normally prefer “‘paved gravel bed’’ to “‘bed with large gravel surface : 
material and gravelly sand subsurface material’ when describing a stream. wom 

_ An example of the strength of the use of * ‘armored gravel bed”” and ‘“‘paved © 

- "gravel bed”’ is a reach of the Yampa River above the Little Snake River in 

a northwestern Colorado. Most the river has an ‘‘armored gravel bed’’ with 5 
_ em-10 cm surface material except in one area where the streambed is “paved” 
by material in the 75 cm range. This larger material was probably transposed 
into the ; Yampa by infrequent high flood flows from the adjacent hills. In this — 
case it is very useful to know that most of the stream is armored, but pavement 
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When describing the ‘material in the stream, ‘the 1 writer — used the term | 
‘armor layer’ to describe the surface layer, and ‘“‘bed material’’ to describe 
the material below the armor layer. In contrast, others have used “bed material” 4 
to > describe a mixture of ‘material in the armor layer and material below the ¥ 


; and subsurfs ace material, is the pref erred v usage e when: describing the bed material. 
The use of ‘‘armored gravel bed’’ and ‘‘paved gravel when referring 
to the channel as a whole and ‘ ‘surface material”’ and **subsurface 
_ material’? when describing the material in the stream bed would go a long way 
in reducing the confusion existing inthe literature. = 
24. Parker, G., and Peterson, A. W., ‘“‘Bar Resistance of Gravel-Bed Streams,”’ Journal 
> om Hydraulics Division, ASCE, Vol. 106, No. HY10, Oct., 1980, pp. 1559-1575. : 


Discussion by Gary Parker,” M. ASCE The 


_ The distinction between armor and pavement offered by the authors senieahien d 
an unfortunate case of inadvertent obfuscation in pursuit of a worthy cause. 

The coarse layer of immobile sediment sometimes produced by os : 
downstream of dams has traditionally been known as armor, not a 


’ The precedents for this are numerous; several are Day (1), Gessler (3), Livesey 
(9), Little and Mayer (20), 4 Ashida and Michive (14), and Garde and Hasan 


: studies. The experiments are characterized by the absence of sediment infeed - 
a recirculation, and thus can only model the phenomenon just downstream | 
aii the context of rivers, at least, the term pavement seems to be of Canadian — 
origin. Kellerhals and Bray (19) have applied the term to coarse surface layers 
7 in streams known to be periodically active. The first person to formally distinguish 
“between the two appears to have been Bray (15), who termed a coarse surface 
a pavement, if surface sizes are fairly common in the substrate, and an armor, | 
if surface sizes are rather rare in the substrate. It is implied, but more or 
_ less obvious from the above, that pavement is periodically mobilized but armor 
“is not. Thus, Bray’s original definition is exactly the opposite of the one offered 
_ The authors may not have been aware that Bray’s original definition has 
taken root in the literature. Gomez (18) specifically associates armor with 
immobility, and identifies pavement as a response to the most recent flood 
_ transport event. ‘Parker (21) and Dhamotharan, Wood, Parker, and Stefan (16) 
_ have gone one step farther. They studied the formation of a coarse pavement 
¥ at low shear stress, but under continuous transport of all available sizes. Both _ 
_ feed and recirculation configurations were used. They showed experimentally ae 
_ that coarse, imbricated surface pavement is so intimately tied up with bed motion 
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DISCUSSION 


_ the coarse pavement ev evolves into an even coarser armor. This line of ing 
a is consistent not only with Bray’s (15) original definition, but also with ox 
‘time- ‘honored usage of the term armor to describe immobile surface layers 
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: an Dam,’ ” Proceedings, 14th Congress, IAHR, Vol. 3, 1971, pp. C30-1-C30-9. err, 
(26. Bray, D. I., ‘‘Generalized Regime-Type Analysis of Alberta Rivers,”’ thesis presented 
; la the University of Alberta, at Edmonton, Alberta, Canada, in 1972, in partial fulfillment eas 
of the requirements for the degree of Doctor of Philosophy. = 
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SAND INSTABILITY Due To Bep Loap 
1¢ following corrections should be made to the original paper: A e 
2038, line 16: Should read ».” instead of “Q 
Page 2045, Eq. 57, right-hand side: Should read Kh, (F’ Kho tanh kh, 
Page 2045, Eq. 60, right-hand side: Should read 


2050, line 20: Should rent Va Deri instead of 


Le 
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Publications, ASCE, 345 East 47th Street, New York, N.Y. 10017. Authors must indicate the 
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to which the paper should be referred. Those who are planning to submit material will expedite 
the review and publication procedures by complying with the following basic ee 2 wid 
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- The manuscript (an original ribbon copy and two duplicate copies) should be euttisambes 
4 on one side of 8-1/2-in. (220-mm) by Il-in. (280-mm) paper. Three copies of all figures and : 


nerally, the maximum length ofa paper is 10,000 word-equivalents. As 
each full manuscript page of text, tables or figures is the equivalent of 300 words. If a particular 
- subject cannot be adequately presented within the 10,000-word limit, the paper should be accompanied — a 
by a rationale for the overlength. This will permit rapid review and approval by the Division 
or Council Publications and Executive Committees and the Society’s Committee on Publications. ; 
Valuable contributions to the Society’s publications are not intended to be discouraged by this 


: ~@ The author’s full name, Society membership grade, and a footnote stating a employment — 
‘must appear on the first page of the paper. Authors need not be Society members. 


gg All mathematics must be typewritten and special symbols must be ‘edinia hiicins The 
letter symbols used should be defined where they first appear, in figures, tables, or text, and 
_ arranged alphabetically in an appendix at the end of the paper titled Appendix. —Notation. i 


6. Standard definitions and symbols should be used. Reference should be made to the lets > 
“published by the American National Standards Institute and to the Authors’ ’ Guide ‘to the Public 
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(220-mm) by I1-in. (280-mm) ‘paper. An explanation of each in the text. 
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